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This thesis is devoted to the development and application of new theoretical 
methods based on the coupled cluster formalism designed for better and more 
general description of linear and nonlinear electro-optic effects in molecules. 

Starting with the general perturbation theory we derive equations for linear 
polarizability and first- and second-order hyperpolarizability. We present both the 
derivation and numerical results for the analytic Equation of Motion Coupled Clus- 
ter (EOM-CC) method, the BW-EOM-CC (Brillouin-Wigner EOM-CC) method 
and the less attractive but more numerically complete energy derivative formula. 

We describe in detail the advantages and disadvantages of including the renor- 
malization terms in the Sum Over States (SOS) expression, widely used for this 
type of calculations, show ways to eliminate some of these terms and discuss their 
numerical importance. 
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Our most powerful method BW-EOM-CC based on the Brillouin-Wigner per- 
turbation expansion is attractive for many reasons. It not only is fully size- 
extensive but also permits calculations of the properties at the frequency of electron 
resonance, where all other methods based on the Rayleigh-Schrodinger perturba- 
tion theory fail. 

We report the numerical results for several interesting systems, including the 
trans-butadiene. Theoretical prediction of nonlinear properties of this molecule 
has become a true test of quality and accuracy of any new method. 
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CHAPTER 1 
INTRODUCTION 

Nonlinear optics is the knowledge about phenomena occurring as a consequence 
of the modification of the optical properties of a chemical system by the presence 
of light. Such phenomena are considered nonlinear because they take place when 
the response of a system depends in a nonlinear (quadratic, cubic, quartic, etc.) 
manner upon the strength of the electric field. 

In 1961 Franken, Hill, Peters and Weinreich [1] demonstrated the first nonlinear 
optical effect, in which ultraviolet radiation at the second harmonic of the ruby 
laser (3471A) was observed when a ruby laser beam at 6942A was focused on a 
quartz crystal. Franken’s experiment confirmed the existence of the nonlinear ef- 
fects and showed that with the optical field strengths available from lasers, effects 
due to the second harmonic generation could become significant. It takes a field 
of a strength of approximately IkV/cm to induce a nonlinear response in a sys- 
tem and no conventional source of light can deliver this intensity. Thus, although 
harmonic generation of electromagnetic waves had been known for a long time, 
physicists needed the power of laser to explore the nonlinear world. Almost at 
the same time, Kaiser and Garrett [2] observed another nonlinear effect produced 
by a laser beam. It was a two-photon absorption giving rise to the second order 
hyperpolarizability . Since then, numerous nonlinear optical phenomena have been 
discovered. They have not only greatly enhanced our knowledge about interaction 
of light and matter, but also created a revolution in optical technology. Nonlinear 
optical properties offer one of the most challenging frontiers for modern quan- 
tum chemical studies. Highly nonlinear materials are being used in many optical 
devices like lasers, optical switches and modulators. There are two aspects of a 
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practical interest in nonlinear optics; The most important of these is frequency 
conversion, in which laser radiation at one frequency is converted, for example 
by harmonic generation or sum-frequency generation, into coherent radiation at 
a new frequency. Because the converted radiation may be at a frequency that 
is not readily available from a laser source, these frequency-conversion techniques 
provide an important means of extending the spectral range of lasers. Converting 
high-energy light from an infrared laser to the visible frequencies would greatly 
increase the density of recorded information on optical media such as CD-ROMs. 
Another aspect is in telecommunication transmission, where information is sent on 
a low energy carrier by placing a dc electric field across the material. This gave 
rise to an extremely fascinating area of optical quantum computing. 

These practical aspects are governed by the interaction of electromagnetic ra- 
diation with matter. 

Optical properties can be expressed in terms of microscopic quantities: po- 
larizabilities and hyperpolarizabilities, which are amenable to quantum chemical 
evaluation. The Hartree-Fock based methods have been known and used for a long 
time (TDHF [3, 4, 5] also known as RPA [6]), but their accuracy is not sufficient to 
be predictive for engineers and chemists. In order to describe hyperpolarizabilities 
more reliably, it is necessary to use sufficiently big basis sets, to introduce electron 
correlation effect, frequency dependence and vibrational corrections. Basis sets 
need the flexibility to describe any polarization of the electronic cloud induced by 
an external electric field, correctly. The effect of electron correlation increases with 
the order of a property; for linear polarizabilities it rarely exceeds 15%, whereas, 
the correlation correction for (3 in the case of ammonia reaches 70% of the total 
value [7]. 

A very important factor is the frequency dependence. The experimental data 
always refer to dynamical processes and therefore relate to frequency dependent 
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properties. The extrapolation to zero frequency often introduces problems for the 
comparison between theory and experiment. 

Vibrational corrections can constitute as much as 10% of the total value of a 
given dynamical property; therefore, some effort is also addressed at including that 
correction [8]. 

Most of the limitations for the electronic problem can be overcome by employing 
correlated methods, together with large basis sets containing polarization and dif- 
fuse functions. So far, there are very few applications of frequency-dependent cor- 
related methods for hyperpolarizabilities: Multi-Configurational Linear Response 
theory (MCLR) [9, 10] and Multi-Configurational Time-Dependent Hartree-Fock 
method (MCTDHF) [11, 12]. Two groups; Rice and Handy [13, 14] and Aiga, 
Sasagane and Itoh [15, 16] derived the formulas for frequency-dependent second- 
order perturbation theory method MBPT(2). For the first hyperpolarizability, the 
last group also extended their quasienergy derivative method for the Brueckner 
coupled-cluster case [17]. Another approach to this problem is the equation of 
motion (EOM) theory [18, 19, 20]. It has been successfully applied to dynami- 
cal polarizability calculations [21], and to the hyperpolarizabilities that occur in 
the Kerr effect by combining analytically computed dynamic polarizabilities with 
finite-field differentiation. 



CHAPTER 2 

MOLECULAR STATIC AND DYNAMIC PROPERTIES 

2.1 Introduction 

The energy of a molecule in an external electric field can be expressed as a 
power series: 

1 1 

E{i) = £'(0) — fiiSi — — 2^ OiijEiCj — — 2^ f^ijk£i£j£k ~ 

i ' ij ' ijk 

^ ^ ^ I ••• (^*^) 

ijkl 

where Ex is the field strength in x direction, etc. and ii, a, j3, 7,... are coefficients 
of the expansion. These coefficients are physical, measureable quantities; dipole 
moment, polarizability, and first and higher hyperpolarizabilities. 

The above equation can be rewritten in terms of the induced dipole moment 
[22] , as the first derivative of the total energy with respect to the external field 
strength i: 



dE{£)\ ^ 1 1 

1 P'i T OLijS-j + ^ . Pijk^j^k T ^ 'JijklSj^k^l + ••• 

i / j 



OE: 






3! 



( 2 . 2 ) 



jkl 



Using this expression we can define the polarizability a as a first-order, linear 
(in external field strength) response property of a molecule in an electric field. The 
first hyperpolarizability /? is then defined as a second-order, nonlinear property, 7 
as a third order property, etc. 

There are two general ways to define molecular properties; 

1 - Energy derivatives, and 

2 - Perturbation expansion. 

The energy derivative method also can be divided in two classes: 
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- Finite field, numerical differentiation, and 

- Analytical derivatives. 

Below we will summarize these methods, paying particular attention to their 
importance in the evaluation of higher-order nonlinear properties. 

Finite field derivatives 

Differentiating Eq. (2.1) with respect to the external field strength, and evalu- 
ating the derivative at the limit ej,€k,ei,€m, ■■■ — >■ 0, we obtain the expression for 
a linear polarizability: 




( 2 . 3 ) 



Similarly for hyperpolarizabilities we get 




( 2 . 4 ) 




( 2 . 5 ) 



and so on. 



Using a simple mathematical formulas for numerical differentiation and ap- 
proximating infinitely small increments of the electric field by finite values we can 
obtain the finite electric field derivative formulas. The dipole moment will be 
then expressed as the finite first-order derivative of the energy with respect to the 
external field strength: 



n(i) = 



( 2 . 6 ) 



C2 - ei 



The best simple method for evaluating numerical derivatives is a central difference 
formula where ci = e and C 2 = — e- This gives 
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Applying the above formula to the linear and nonlinear polarizabilities we obtain 



Oiij — 



_ ~ + Ei(—2ej) — 2Ei(0) 



2e. 



4e] 



(2.8) 



f^ijk 



^ij{ ^k) 



26 . 



^ijk{^l) Pijki, 

im - ^7 



(2.9) 



( 2 . 10 ) 



The numerical differences method is a very convenient method as it delivers 
values of the properties at any order at which the single point energy can be 
calculated. It can be thus extremely useful serving as a good benchmark and 
assuring accuracy of a new method. 

Drawbacks of the finite field method are on the other hand numerous. First of 
all it is very expensive-a single tensor element of the second hyperpolarizability 
7 requires 5 separate calculations at different electric field strengths. Also, the 
choice of the field strength is a trade-off between the accuracy of the result (the 
smaller the field is, the closer to the full derivative the result is) and the magnitude 
of an effect (large field causes a bigger response and more significant digits in the 
answer). Also, it is very difficult if not impossible to calculate any frequency 
dependent properties, as it requires an analytical frequency dependent code for 
lower-order derivatives. 

Derivatives as perturbative energy corrections 
Let us define a perturbed Hamiltionian as 



H{e) ^ Ho + eV, 



( 2 . 11 ) 



where Hq is the unperturbed Hamiltonian, and V is the perturbation. For the 
properties of interest, F is a electric dipole moment in x,y or z direction. The 
corresponding Schrbdinger equation with e-dependent wavefunction, Hamiltionian 
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and eigenenergy can be written as 



H{e)tp{e) = E{s)xp{e). 



( 2 . 12 ) 



Differentiating the resulting expectation value of H{s) with respect to the field 
strength, we get 



E{£){ip{e)\ip{e)) = 



(2.13) 









ds 



de 



mme)) 



dH{e) dE{e) 



de 



de 



IV’(£» 



(2.14) 



Two first terms of the above equation vanish, since as e — )• 0, H{e) — )• Hq, E{e) — > 
Eq and we have {Hq — Eq) \ipo) — {fpo\{Eo — Eq) = 0. This leads to 



dE{e) 

de 






('0o| ?T^l^o)|e =0 — (V'ol^lV’o) = 



de 



(2.15) 



Taking the second derivative of the expectation value (2.14) gives (assuming the 
wavefunction to be real) 



d^Eje) 

de'^ 



£=0 



= 2(^, 



0 



= 2<V>o|V'|!A‘'>) = 2Eg>p^. (2.16) 



And in general, 
d^^^E{e) 



de 



n 



de ' de(^-d ^ 



{'ipo\V\ip''^ = (^O'^HSPr 



(2.17) 



Hence the perturbation theory may be straightforwardly used in evaluation of 



molecular response properties, like (hyper)polarizabilities. 



2.1.1 First Hyperpolarizability, P 

The first polarizability /3 corresponds to a second-order, nonlinear process, and 
it can be categorized by frequency, intensity and phases of the field components. 
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The second-order nonlinear tensor element is expressed as 



(2-18) 

where u >3 — uJi i: LO 2 , and i,j,k are the Cartesian components of the interacting 
fields and polarization waves. In general, each component is associated with a 
frequency component in the argument of (3. Assuming the Kleinmann symmetry 
[23], i.e. choosing the incident frequencies far from resonance, Cartesian indices 
and frequency components can be permuted together. There are two physical in- 
terpretations of /3; in one of them, the field associated with the first beam with 
frequency Ui alters the refractive index of the medium and the propagation char- 
acteristics of the second beam with frequency U 2 is modified accordingly, is phase 
modulated and exhibits sidebands at frequencies u)\ + ui 2 and ui\ — uj 2 - If both in- 
coming beams have the same frequency, a harmonic overtone at 2uj is created. The 
other interpretation describes the nonlinear effects as resulting from nonlinearities 
in the polarization response to incident fields at various frequencies. A medium 
is composed of a set of classical anharmonic oscillators, describing, physically, an 
electron bound to a core. Here, new frequency components of the polarization at 
ooi ± cd 2 , 2uJi, 2 co 2 and 0 appear through quadratic interaction of the field with 
oscillators via the anharmonic terms. 

In the case where cui = u>2 = lv, we have one of the best known effects in non- 
linear optics; second harmonic generation, SHG expressed as P{—2u-,u,uj). In this 
process, the energy is redistributed between the fields as a result of interaction of 
the waves with the medium. To have efficient redistribution of energy, momentum 
must be conserved. This condition is known as phase matching. The implication 
here is that energy can only be transferred from the fundamental to the harmonic 
field, as the waves propagate through the medium, as long as they maintain the 
required phase relationship. 
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If Wi = 0 and —003 = C 02 — , vfe speak about an electrooptic effect (EOE), 

where a dc static field is applied to the medium, altering the refractive index, and, 
hence, the propagation characteristics rather than the frequency of tight. This is 
written as 

(3ijk{-u};0,uj) (2.19) 



The last second-order nonlinear effect is an optical rectification (OR), 

f3ijk{0]-u,u) (2.20) 

for which tJs = 0 and coi = —u> 2 . Here, two beams with frequency u interact 
to produce a dc polarization of the medium. Assuming validity of the Kleinman 
symmetry [23], we have 



/^a6c( CJ,0,Ct^) Pcabi,^) ^?^)? 



( 2 . 21 ) 



which shows that EOE and OR are closely related. 

From Eq. (2.17) we note that to calculate /?, one needs to know the second- 
order perturbed wavefunction. Recalling the RSPT formulas; 



^(2/c+i) ^ ^ 



( 2 . 22 ) 



where 



k — l 



^/;W) = i?o (R - £:O|.0(fc-O) 

1=2 



(2.23) 



Here 

and 



is the n-th order wavefunction, Fis an external perturbation potential, 
Rq is the resolvent operator; 



flo = E IV'.' 



(0) 






(0) 



(2.24) 



10 



The first hyperpolarizability can be viewed as 

(3 = (2.25) 

Using Eq. (2.23) the above expression can be rewritten as 

dijk = 6(4°)|Ki?o (H - (H)) i?oK|4°’) ■ (2.26) 

The full expression for a specific tensor element of the second harmonic generation 
then becomes [24] 

2cj,cj,cj) 2cj, u^) “f” Eijf^(^^co ^ Lo^ T 2u^, u)'^ T 

Fikj{2u,uj) + Fkij{—2uj, —uj) + Fkij{2uj,uj) (2.27) 



with the frequency dependent moment defined as 



Fijk{nu,muj) = ^ 

pq^O 



(Fp + nuj) {Fg + mu) 



(2.28) 



{Vx) being the expectation value of the perturbation, (V'o| 14 |^o) • 



2 . 1.2 Second Hyperpolarizability, 7 

In general, third-order nonlinear properties can be viewed as resulting from 
four-wave mixing. There are numerous third-order processes that can be distin- 
guished on the basis of frequencies of the output and input beams. The most 
important case is a Third-Harmonic Generation (THG), where all three incoming 
waves have the same frequency u, and the resulting response is a beam at the 
frequency 3a;; 

7q/fe( (-3a;;a;,a;,a;) . (2.29) 



If one of the frequencies has a negative phase, we observe an effect called the 
Intensity-Dependent Refractive Index (IDRI), and the process is described by the 
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following expression: 



lijkl {-UJ-UJ,UJ, -Oj) . 



(2.30) 



Incident frequencies do not necessarily have to be the same; an example is one of 
the most important third-order effects, called the Optical Kerr Effect (OKE): 



lijkl <^1,^2-, —^ 2 ) ■ (2-31) 

Introducing a static dc component of the incident field, we can observe 4 more pro- 
cesses: direct current-induced SHG (dc-SHG), Electro-Optic Kerr Effect (EOKE) 
and direct current-induced Optical Rectification (dc-OR). 

From Eq. (2.17) we have 

7 = 24E;(^’ = 24(V’'“^|K|V'(^)) = 24(^(^)|K - 

(2.32) 



Inserting Eq. (2.23) into the above expression we obtain 



jijki = 24{'iP^^^ViRo {Vj - E^'^) Ro {Vk - Ei'^) 
24£;(2)(v,(«)|I4R2^|^(o)). 



(2.33) 

The second term is a normalization part cancelling unlinked diagrams arising in 
the former term. If we include only linked terms in the above equation, we can 
leave the renormalization term out: 



jijki = 24(7AC’)|KRo {Vj - Ro {Vk 



Ei^^)RoVi\7p^^^)L. 



2.2 Description Using Hartree-Fock Methods 



(2.34) 



Theoretical description of linear and nonlinear processes in molecules using 
analytical (not finite-field derivatives) methods began almost simultaneously with 
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the development of time dependent field wavefunction. The fully analytical formu- 
lation of the ab initio frequency dependent polarizability and hyperpolarizabilities 
dates back to 1986, when H. Sekino and R.J. Bartlett [3] derived and evaluated 
these properties using time dependent Hartree-Fock theory (TDHF). They utilized 
the recursiveness of the TDHF theory and developed formulas that were applicable 
to any order. This method also incorporates the 2n + 1 rule for energy expres- 
sion [25], stating that the 2n + 1 energy order can be computed by using just the 
nth-order TDHF wavefunction, which can offer an additional simplification. 

Starting from the Schrodinger equation in the time dependent formalism. 



FC - = SCe , 

dt 



(2.35) 



and expanding the time dependent coefficient matrix C by orders with respect 
to the external oscillating perturbation and limiting the derivation to a single 
monochromatic and static fields Sekino and Bartlett obtained the final expressions 
in terms of density matrices: 



«(-o;;a;) = -Tr {h<^'D(±u;)} , 

P{-2uj;uj,uj) = -2Tr , 

7 (— 3a;; w, a;, a;, w) = — 6Tr ±o;, ±c<;)| ... (2.36) 

D is the density matrix defined for first-order as 



D(±a;) = C(±a;)nC(°'t + C^°^nC{Ttoy , (2.37) 



and for the second-order as 



D(±a;, ±io) = C(±a;, ±o;)nC(°^t + C(±a;)nC(Ta;)^ + C(°)nC(Ta;, . (2.38) 

In general, for any meaningful order n, the TDHF hyperpolarizability can be ex- 
pressed in terms of the n-th derivative of the density matrix with respect to the 
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electric field strength E: 

= -Tr|HC)^D| 

= -n!Tr {h(^)D(±o;,...)} . 

The static (time independent) formulation of the coupled-perturbed Hartree-Fock 
(CPHF) obeying the 2n + 1 rule had been implemented by Dykstra and Jasien [26] 
a few years earlier in 1984, but it is the frequency or time dependent processes that 
are of the highest importance in chemistry and physics. Thus we will concentrate 
on methods that include the frequency dependence. 

Unfortunately, the accuracy of the TDHF based methods is not sufficient to 
be predictive for engineers and chemists. In order to describe hyperpolarizabilities 
more reliably, it is necessary to use sufficiently large basis sets and introduce elec- 
tron correlation effect. Basis sets need the flexibility to describe any polarization 
of the electronic cloud induced by an external electric field, correctly. The effect of 
electron correlation increases with the order of a property; for linear polarizabili- 
ties it rarely exceeds 15%, whereas, the correlation correction for in the case of 
ammonia reaches 70% of the total value [7]. 

2.3 Description Using Correlated Methods 
2.3.1 Many-Body Perturbation Theory Based Methods 

The first known correlated, single-reference, fully analytical treatment of fre- 
quency dependent hyperpolarizabilities dates back to 1991, when Rice and Handy 
[13] formulated and implemented the second-order Mpller-Plesset perturbation the- 
ory expressions. To achieve their goal, the authors used a concept of pseudo-energy 
W{t) defined to be the time-dependent eigenvalue of the operator H — i[d{dt)\. 

f H - $ = W(t)$ 



(2.39) 
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Here, solution of this equation is the same as solving the standard time- 
dependent Schrodinger equation 




(2.40) 



where 






(2.41) 



W{t) was then expanded in a power series in E{lo) = E{<d)-\-\/2E{uj) 
with W (0) being the static eigenvalue. The partitioning of the Hamiltonian is then 
applied, where the unperturbed part is the time independent Hamiltonian Hq and 
the perturbation H' is defined as —i[d/{dt)\. 



Ho = E-i 



dt 



(2.42) 



H' = H - E 



(2.43) 



where 

F = Fo- (2.44) 

From here, the first derivative of the quasi-energy W with respect to the external 
field strength (time-dependent dipole moment), second derivative (polarizability) 
and third derivative (first hyperpolarizability (3 is evaluated as a complicated and 
not recursive expression [13, 14]) as corrections to the unperturbed (static) values: 



d d d'it 

aij{—uj-,u) cos{uJt) = 0 (^( 0 ; 0) cos(a;t) + *^ 77 , 



dEddt' 'dE. 



Ul] 



0^ d 

(3ijk{-OJ\0,Uj) COs{ujt) = /?jjVt(0;0,0)cOS(o;t) -I- ) 



dEidEj dt dEi^k 



(2.45) 

(2.46) 



The numerical results presented by Rice and Handy on ammonia and formalde- 
hyde prove that it is desirable to determine the frequency dependent contribution 
to the (hyper)polarizability at the correlated rather than SCF level, especially in 
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cases, where the frequency dependent contribution is significant, such as the SHG 
process. 



2.3.2 Linear and Quadratic Response Methods 

Another way to describe the linear and nonlinear properties is the response 
function [27], or propagator formalism [28, 29]. The response method can be 
applied to any wavefunction, and therefore there have been many different for- 
mulations of the method. Response theory methods for calculations of molecular 
properties have recently been vastly improved. These methods are based on both 
non-correlated SCF and correlated and multi-configurational SCF (MCSCF) ref- 
erence wave functions. 

The response functions are defined by considering the expansion coefficients of 
the expectation value of a (time-independent) perturbation operator. In case of 
optical properties of interest, this operator is the dipole moment operator // [30]. 



= (T(t)|/r|T(t)) = (^l/i]T) -h J {{ijL]V{u;i)))^, exp{-iuit)dui 

+ 1/2 J J {{fJ--,V{uJi);V{u2)))u;i,u2 exp(-io;i - iu2)tdujiduj2 

+1/6 J J j {{p.\V{iOi)\V{uj2)\V{uJz)))u,u>^2,u^i^w{~i^i-i^2-i^z)td(jOiduj2duJz . 

(2.47) 



Comparing Eqns. (2.2) and (2.47) one can readily see that the response funct- 
tions ((/x; V (a;)...)) define the properties of interest. The first, linear response func- 
tion defines polarizability a and can be evaluated as a Sum-Over-States expression 
1301 : 






LJl 



u 

k 



(To]V;(a;i)|T("))(v[rWl/Xi]To) 



}• 






^1 + ^k 



(2.48) 
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Similarily, the quadratic response function defines the first hyperpolarizability (5: 



(o;i + <jJ2 — LOk) (<^2 ~ <^n) 

(^o|yfc(g;2)|^^"^) -4n(^o|V^,(^i)|^o)) 

(uJi + (jJ2 + iOk) (^2 + ^n) 

(wi + Uk) (u^2 — UJn) 

(2.49) 



)) (^(")|14(g;2)|^o) 



E{ 

k,n 







assuming that the summation indices run over both ground and excited states. If 
we exclude the ground state from the summation, the resulting expressions contain 
several more terms, called ’’renormalization terms” [31, 32]. 

SCF linear and quadratic response functions were derived for the first time 
in 1961 [33], while the first multiconfigurational (MCSCF) formulation of linear 
response was presented by Yeager and Jprgensen [34] and independently by Dal- 
gaard [35]. A more general formulation of the MCSCF linear, quadratic and cu- 
bic response functions was given by Olsen et al [29]. Response functions have 
also been derived for other wavefunctions: linear response has been obtained for 
the MBPT(2) wavefunction in the polarizator propagator formalism [36]. Also, 
a coupled-cluster wavefunction has been used extensively to determine response 
functions. The linear coupled-cluster response equations have been formulated 
by Monkhorst [37] and Dalgaard and Monkhorst [38] and extended to quadratic 
response and applied to calculations of polarizability and hyperpolarizability by 
Koch and Jprgensen [39, 40]. 



CHAPTER 3 

COMPUTATIONAL METHODS USING COUPLED-CLUSTER 

WAVEFUNCTION 

3.1 Rayleigh-Schrodinger Perturbation Theory Based Methods 

In this chapter we will derive the expressions for linear polarizability and first 
and second-order hyperpolarizabilities in the presence of a static electric field. 
Consequently we will present the derivation of the time or frequency dependent 
first polarizability and show that the only difference between static and dynamic 
expressions is the presence of the frequency factor in the resolvent. The subsequent 
treatment of higher-order properties and new methods will be based on the time- 
independent formalism with frequency components inserted directly into the static 
expressions. 

We will start by defining our Schrodinger equation as 

{Ho + V) '^ = E'i . (3.1) 

where Hq is the unperturbed Hamiltonian satisfying the equation 

= iEo$o (3-2) 

and V is the perturbation (response of the system to the external electric field at 
frequency a; = 0). 

To define our methods we will rearrange the equation (3.1) by adding an arbi- 
trary factor A: 

(/fo + U + A) ^ = (^ - A) (3.3) 
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Considering Eq. (3.2) and assuming intermediate normalization ((<l>ol'I') = 0) 
we may write Eq. (3.3) as 



(A - //o) ^ = (E - E - A) ^ (3.4) 

Assigning different values to the parameter A leads to various types of a pertur- 
bation expansion. However, we have to restrict the choice of A to values different 
from the eigenvalues of Hq in the ^ — $o space (also known as the (^-space, or- 
thogonal to the subspac^ P which we choose to be ^'o)- The reason for such a 

/N /S 

restriction is to allow for the inverse of Q (A — Hq) Q'if in a modified Schrodinger 
equation 

Q{\-Ho)Q^ = Q{V-E + X)<i> (3.5) 



obtained from Eq. 
commutes with Hq. 
operator [41] 



/N /s 

(3.4) by using the indempotency of Q and the fact that Q 

/>. /N 

The inverse of Q (A — Hq) is also defined as the resolvent 



Ro{x) 



Q 

X-Ho 



r Ei) 




Substituting ^ = $o + <3^ using the above definition of the resolvent we can 
rewrite Eq.(3.5) to read: 



= $0 + Po(X) (V - B + X)^ 

= $0 + Po(X) (V-B + X) + Po(X) (V-E + X) = --- 

OO 

= ^{/?o(A)(E-^ + A)r^o (3.7) 

1=0 

If we define the total perturbed energy as £” = E'o -t- AE, we can then evaluate 
the perturbation energy as 

OO 

AE = J^(^oiV{Po(X) (V-E-h A)}'|$o) • 

i=0 



(3.8) 
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The Eq.(3.8) allows us to calculate the perturbed energies knowing only the 
ground state wavefunction. 

All the methods utilizing the perturbation expansion to treat excitations and 
properties were based on the Rayleigh-Schrodinger perturbation theory (RSPT), 
for which we define the lambda parameter (3.4) to be the ground state energy 

A = £’o- 

To apply RSPT to our purpose we will use the order-by-order expansion, fol- 
lowed by the zeroth-order expansion of resulting expressions, in order to get equa- 
tions that are explicitly dependent only on the ground state wavefunction. We 
start with expanding the perturbed wavefunction and energy in a power series 
around the unperturbed values: 

^ = ^ + x = + . . . , 

E = + AE = + . . . . (3.9) 

By using the coupled-cluster parametrization for the ground and excited states 
(EOM-CC states), we can redefine the wavefunction 

^cc = e^|0) (3.10) 

and our Hamiltonian 

H = e~'^HNe^ =Hm + kV . (3.11) 

We will refer to .5^ as a similarity-transformed Hamiltonian. is the normal- 
order Hamiltonian (/fjv = H — (0|i/|0)). The ground state |0) is a reference 
determinant. The similarity transformation is non-unitary and therefore the re- 
sulting Hamiltonian is no longer hermitian. Hence, the bra and ket wavefunctions 
are no longer simply orthogonal or each others conjugate, but instead they form a 
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biorthogonal basis [20]: 



{tpo\LkRi\'4’o) = C6ki 



(3.12) 



In the EOM method, we choose Lq = (1 + A), Rq = 1 and C=1 for the ground 
state, where A is a deexcitation operator complementary to the excitation operator 
T [42], We also partition our space into three subspaces: the ground state (Fermi 
vacuum) |0)(0|, the space of single and double excitations |h)(h| (because we work 
within the CCSD approximation; had we used a more rigorous CCSDT method, 
this subspace would contain triple excitations) and the remainder - all higher 
excitations |f)(f] Together all three subspaces form a full space spanned by the 
operator H. 

It is noteworthy to mention that the above formulation satisfies the generalized 
Hellmann-Feynman theorem [42]: 

^{0\{1 + A)H + Kt/|0) = (0|(1 + A)F|0) (3.13) 

QjKj 

Rewritting the Schrbdinger equation (3.1) with 



^ ys 



H = Ho + kV, 



(3.14) 



we get 



(Ho + kV - + Av2<J,(2) + ^ = 0. 



(3.15) 



Equating coefficients of equal powers of k gives 



m—2 






1=0 



m—1 



^^( 0 ) _ fjA ^ E^rn-l)^(l) 



1=0 



(3.16) 
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for k"" (mth-order equation). After applying the ground state left eigenvector and 
expressing the perturbed wavefunctions in terms of the zero order wavefunction 
we obtain 

($o| (^0 - = (<E>o|b^|«J>o) , (3.17) 

($o| {Ho - = ($o|iy|$(i)) + £( 2 )|$( 0 )^ 

= {^o\WRo{Eo)V\^o) + , (3.18) 

(<l>ol {Ho - E^^^) |$(^)) = + E^^^ + 

= {^o\WRo{Eo)WRo{Eo)V\^o) + E^^^ + E^‘^^^o\Ro{Eo)V\^o) , 

(3.19) 

($o| {Ho - + E^^^ + 

= {^o\WRo{Eo)WRo{Eo)WRo{Eo)V\^o) + E^^^ + E^^\^o\R{E)oV\^o) 

+E^^H%\Ro{Eo)VRo{Eo)V\^o) + E^^^^o\VRo{Eo)Ro{Eo)V\^o) ■ (3.20) 

In the above equations W = V - E^^"> and Ro{Eo) is the RS resolvent operator 

Ro(Eo) = E mmEo - (3.21) 

u 

defined for the non-diagonal H Hamiltonian. The RSPT resolvent contains only 
the zeroth-order energy and thus evaluating equations (3.17) to (3.20) is straight- 
forward and does not involve any iterating procedures. 

The coupled-cluster perturbed wavefunction |‘P(n)) can be written in terms of 
the unperturbed ground state wavefunction with the excitation operators being 
actually perturbed: 

^ 7 ^( 2 ) ^ 7 ^( 1 ) ^ 




1 ^( 4 )^^^ _ 2 ^( 4 ) _|_ 2"(3)2'(1) _|_ 2 '( 2 ) 2 "{ 2 ) _|_ rj,(2)rjn{l)rjl{l) 2’(^)J'(1)J'(1)X'(1)|$q\ 

(3.22) 



Applying the above results to the zero order expansions in Eqs. (3.17) to (3.20) 
and noting that (<ho| [Hq — T(”^|d>o) = 0, ($o| = (0|(1 + A) and |<I>o) = |0) 



we get 



= (0|(1 + A)1/|0) , 



(3.23) 



= (0|(1 + A)WRo{Eo)V\0) + (0|(1 + A) (^o - E^^^) , (3.24) 



(0|(1 + A)WRo{Eo)WRo{Eo)V\0) + 

E(^)(0|(1 + A)i?o(Eo)E|0) + (0|(1 + A) (^0 - E^^^) + 

(0|(1 + A) [Ho - E^°^) , (3.25) 



= (01(1 + A)WRo(Eo)WRo(Eo)WRo(Eo)VlO) + 

E("^(0|(1 + A)i?o(Eo)E|0) + E(2)(0|(1 + A)Ro(Eo)lVRo(Eo)VlO) + 
E(^^(0j(l + A)WRo(Eo)Ro(Eo)VlO) + 

(0|(1 + A) (Ho - E^^^) T^^)t^^)t^^)t^'^)\o) + 

(0|(1 + A) (^0 - + 

(0|(1 + A) (Fo - E^^^) + 

(0|(1 + A) (^0 - T(2)t(2)|0) . 



(3.26) 



The above equations contain three classes of terms; the leading term is a 
sum-over-states-like term which is the most significant numerically and has its di- 
rect equivalent in the polarizator propagator formalism [27, 28, 29]. The second 
category of terms is the renormalization terms coming from the perturbation ex- 
pansion. They contain unlinked diagrams of type * (0|/(7?, E)|0) which cancel 
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unlinked terms showing up in the leading term, due to the structure of both the 
perturbation operator {V - and the Cl-like operator A. The last part of the 
renormalization terms is the ”multiple-T”, higher-order terms coming from the 
expansion of the perturbed coupled-cluster wavefunction. They describe effects of 
higher (triple, quadruple etc.) excitations that occur in the leading SOS-like term 
but are not adequately accounted for in the CCSD partitioning of the full space. 
The numerical significance of the renormalization terms is discussed elsewhere in 
this dissertation. 

We can now utilize Eqs. (2.25) and (2.32) and write down the RSPT equations 
for linear polarizability as well as for the first and second hyperpolarizability for a 
case where the higher-order renormalization terms are neglected: 

a{0; 0) = 2(0|(1 + A)Wi?o(£^o)R|0) , (3.27) 



/?(0;,0,0) = 6{0\{l + A)WRo{Eo)WRoiEo)V\0) + 

a(0;0)(0|(l + A)Ro(^o)R|0) , 



(3.28) 



7(0; 0, 0, 0) = 24(0|(1 + A)WRo{Eo)WRo{Eo)WRo{Eo)V\0) + 

/?(0; 0, 0)(0|(1 + A)i?o(^o)P|0) + a(0; 0)(0|(1 + A)Ro{Eo)WRo{Eo)V\0) + 
a(0;0)(0|(l + A)Wi?o(^o)Ro(£^o)^|0) . 

(3.29) 



3.1.1 Time-Dependent Derivation of the First Polarizability Expressions 



Properties can usually be evaluated from the derivative of energy associated 
with a wave function, subject to an external perturbation. For time-independent 
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perturbations this is straightforward. For example, dipole-moments, static po- 
larizabilities and hyperpolarizabilities are defined as the first, second and third 
derivatives of energy which are proportional to the zeroth-, first- and second-order 
responses of a molecule to a static electric field. However, a straightforward ex- 
tension of energy derivative methods to time-dependent perturbations is not as 
apparent, because the corresponding energy is not uniquely defined. Hence, the 
starting point of our derivation is the time dependent dipole moment which is well 
defined and unique. 

The time-dependent dipole moment is given by 



nit) = (4^ (t) (t)) 



(3.30) 



where |'F (t)) satisfies the time-dependent Schrodinger equation 






(3.31) 



and is assumed to be normalized at time t = 0, and hence at all times. To 
define polarizabilities, hyperpolarizibilities and so forth, let us consider the field 
dependent Hamiltonian 



H{t) = i/o + XI (<^a) /^ae 



— iuJat 



a 






(3.32) 



where //q is the usual molecular Hamiltonian, fia denotes the three Cartesian com- 
ponents of the dipole operator, while Sa {uJa) indicates the respective components 
of the frequency-dependent field strength. (Every integration, either in time or 
in frequency is accompanied by a factor . We absorb this factor in our 

integral sign to facilitate the notation. 

At time t = 0 we assume that the state (0)) is the ground state (or possi- 
bly another stationary state) of the field independent molecular Hamiltonian i/o- 
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Under these conditions the polarizability and higher polarizabilities (with respect 
to the state at time t = 0) can be defined as the expansion coefficients of the 
time-dependent dipole moment 



11 A (t) = ^J‘A^ + Y^ / duJaO-Aa £ («^a) 6 

a 



— iuJat 



-t- 



^ r 

— ^ y dUadUbpAab {^a, ^b) £ (^a) £ (^^6) -f- ... 

a, 6 



(3.33) 



Here, we adopt the convention [15], using upper case letters for the component 
under consideration and lower case letters to indicate components that are summed 
over. For future reference, the Fourier transform of /ia (t) reads 



Ha (o;) = 



‘OO 



Ha (t) 



— OO 



H^^5 (a;) -h ^ / dbOaOCAa M £a (t^a) 5 {uj - Ua) + 

a 

/ OO 

doJg^duJbPAab ^b) £a i^a) £b (^b) d ^b) d" ... 

a,o 



(3.34) 



It follows that the polarizability can be obtained as 



OiAB 



f°° Oha {uj) 



E=0 



(3.35) 



where E indicates all components of the field strength. The first hyperpolariz- 
ability is given by 



/ OO 
-OO 



du) 



8ha (t^) 



Oeb {^b) dsc i^c) 



E=0 



(3.36) 



The wave-function (t)) that defines the time- or frequency dependent dipole 
moment satisifies the time-dependent Schrodinger equation corresponding to the 
complete, field dependent, Hamiltonian. This wave function (and the time-dependent 
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dipole moment) can be expanded in orders of the field, and this will yield expres- 
sions for the polarizability and hyperpolarizabilities. If we assume that all eigen- 
states of the molecular Hamiltonian are known (at least formally), we can obtain 
in this way the familiar sum over states expressions for the polarizability and so 
forth. Alternatively, we can assume a parameterization for the complete field- 
dependent wave function, which is then subject to a perturbation expansion. This 
leads directly to equations for the perturbed amplitudes, which in turn define the 
quantities of interest. One possible parameterization is, of course, the exponential 
coupled-cluster ansatz, which is the focus of this thesis. The above approach to 
polarizabilities (or hyperpolarizabilities and so forth) does not apply directly to 
the coupled-cluster formalism, however. The reason is that the time-dependent 
wave-function in the CC formalism is not normalized to unity but rather satisfies 
intermediate normalization 



^cc (t)) — — >■ ('^‘ol^cc (t)) = 1 



( 3 . 37 ) 



The exact CC wave-function is proportional to (t)), where the proportionality 
factor is in general complex. Let us write 



= \^cc (i))e-*'"« 



( 3 . 38 ) 



where the normalization factor is defined as 









( 3 . 39 ) 



Substituting this expression in the time-dependent Schrodinger equation, and mul- 
tiplying by we obtain 



it)) + it)) = Hl'^cc it)) 



dt 



dt 



( 3 . 40 ) 
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If we multiply in addition by e we find 



"?^|ci> + ^|o) = «W |o) 



dt 



dt 



(3.41) 



where H (t) = e It may be appreciated that the equations for T are 

completely decoupled from the equation for the phase 



{q\H {t) - i^\Q) = 0 



dt 



dip 

dt 



= (0|H(t)|0) 



(3.42) 



The analogy with time-independent CC theory is evident. The quantity ^ is 



called the quasi-energy, and it has been observed that polarizabilities and so forth 
can be obtained from differentiating this quantity [38, 37, 15, 13, 16]. The relation 
between the quasi-energy and the quantities of interest is not very transparent in 
our opinion; however, so we prefer to not use the quasi-energy. The time-dependent 
dipole moment in terms of the CC wave-function is defined as the expectation value. 






ei'Pit) {<l! CC (t) (t))e-*VW 



= (^cc (0 IAaI'I'cc {t)) 



(3.43) 



The state 



{^cc (t) I = 



(^cc (i) 



(^cc (t) \^cc {t)) 



(3.44) 



by definition satisfies the normalization condition 



{^cc (t) l^'cc {i)) = 1 



(3.45) 



and it can be conveniently parameterized as 



i^cc it) I = ($o| (1 + A it)) (3.46) 

where A (t) is a de-excitation operator introduced in analogy with CC derivative 
theory [43, 42). The normalization condition is satisfied trivially for all values of 
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the parameters determining the operators A (t) and T (t). The state (^'cc (0 



is related to the true Schrodinger bra as (t) \ = {'i! cc {t) \ the time 

evolution of this state is of course determined by the Schrodinger equation 






(3.47) 



If we substitute the parameterization for (T (t) \ in the Schrodinger equation and 
multiply by and from the right, we find 



(0|(1 + A(t)) 



H{t) 



^ 

* dt dt 



+ i(0| 



dA jt) 
dt 



0 



Projecting on |0) and \q) = g|0), we obtain 



d(f 

dt 



= (0|5(()|0) + (0|A(() 
= (0|ff(f)|0) 



H(t) 



dt 



| 0 ) 



(3.48) 



(3.49) 



and 



(0|(1 + A(t)) [H {t),q\ 10) +i(0|^|9) + 
(0|(1 + A(t)) J^(t)-^^-^||0)=0. 



(3.50) 



The last term in this equation vanishes because the CC equations are satisfied 
(this is true also if T is truncated). Let us summarize our basic equations obtained 
so far: 



/r (t) = (0| (1 + A (t)) e = (0| (1 + A (t)) jl {t) |0) , (3.51) 



(9l li/W-i^||0) = 0, 



(3.52) 



(0| (1 + A(t)) [h {t),q] jO) + i(0|^|9) = 0 . 



(3.53) 
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In addition, we have an equation for the phase factor that establishes the 
connection between our CC states and the true (normalized) wave functions 



dip 



= {0\H{t)\0) 



( 3 . 54 ) 



However, this quantity will not play any further role in our derivation. To obtain 
a perturbation expansion for the time-dependent dipole moment, we assume a 
perturbation expansion for the cluster operator T and the operator A. 



OO 



T{K,t) = {t) 



n=0 



OO 



A(K,t) = (t) 



n=0 



To keep the notation compact we also introduce an expansion for 



OO 



H{t) = Y (t) 



n=0 



( 3 . 55 ) 



( 3 . 56 ) 



( 3 . 57 ) 



The components are defined in terms of Hq,V and the various for 

example 



^(0) ^ 



^(1) ^ 






ff{n) _ fj{0) J^in) jj{n/n-l) 



( 3 . 58 ) 



For future convenience we have introduced to denote all terms in that 

can be obtained by connecting cluster operators of order < m to H = Hq + V . 
This notation allows us to write explicit equations for for n > 0 



{q\ (t) 



QJ'in) 

dt 



|0) = -(9|^("/”-i)(t)|o) 



( 3 . 59 ) 
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These equations can be solved hierarchically. Similarly the equations for can 



be written 



{0|A("> it) [ff<"'.9]|0) + i{0| — 1,) = 



dt 



-(0| (l + |0) - ^(0|Af^)(t) |0) . (3.60) 

k=l 

Again these equations can be solved hierarchically. The equations for A^"\ require 
the solution of perturbed amplitudes T up to and including The zeroth order 
equations are the usual time-independent T and A equations 



= 0 



(3.61) 



and 



(0|A[.^W,g]|0) = -(0|^W|9) 



(3.62) 



The higher order equations are best solved in the frequency domain. If we introduce 
Fourier transforms 



T{uj)= I (t) dt 



(3.63) 



the equation for (a;) can be written 



{q\ (o;)] (lo) |0) = (a;) |0) . 



(3.64) 



With the operator (^1 |0) + lodpq we associate its inverse, the resolvent 

/s 

BP (u) , that only acts in the space of excited determinants. Hence we can write 



{q\T^^^ (u) |0) = -{q\R^ (co) (w) |0) . (3.65) 

The resolvent has the very useful property that when acting on a connected 
operator in its right, the result is again connected. We will see that the structure 
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of the whole approach will be intrinsically connected. The quantity 



|o) = f (t) \ 0 )dt (3.66) 

is in general a complicated expression that involves the Fourier transform of a prod- 
uct of functions of time (T-amplitudes, a time-dependent perturbation etc.). Each 
of these elementary functions have a Fourier transform, and the total transform 
can be written as a generalized convolution, e.g. 



/ {t) g (t) h (t) = 



j j j dujiduj2dujz J f {ui) g (002) h (cos) 

j J J dxJidoj2di03f (uji) g {U2) h (uJs) 6 {uj - ui - U2 - UJ3) . (3.67) 



We note that these type of expressions have precisely the same form as the 
polarizibilities and so forth, and this will therefore provide an easy connection. 
The equation for takes the form 



(0|A<“'H[5<“l,9]|0)-a,(0|A(">(c.;)|,) 
(0| (1 + A<»>) [/?<“) H,g]|0) 



n—1 



- $:(0|A('=)(a;)[^("-^)(a;),(7]|0). 

k=l 



(3.68) 



Using the resolvent, the solution can be written 



n—1 



{0|A>”> (u,) Ip) = - y ^{0| + A<‘) (a.)) (a,) , 5 ] |0)/?;,(-u-) . (3.69) 



9 fc =0 

The above equations are not the most convenient for application. One may 
manipulate further starting from the expression for the frequency dependent dipole 
moment and eliminating the frequency dependent A terms. 
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As mentioned before, the dynamic polarizibility is related to the first-order 



term in the dipole moment, 



oiAa (w) = 



-h 



= (0| (l + A<">) |0> 

(0|AW(^)A“’|0) 



(3.70) 



where 



{,|t('>(u-)|o) = -(?|h„Hp<")|o) 



(3.71) 



and 



(0|Al‘) (u) Ip) = (0| (l + A(»)) [5<') (cp) , ,] |0 >k; (-0-) 



(3.72) 



Using 






r.O |n\ T^(l) 



(3.73) 



and 






(3.74) 



we find 



aA. M = (0| (l + A(»>) [p?>, ri'l (a;)] |0) + (0| (l + A(»>) [p<“>.T<'’ |0> 



-f 



(0| (l + A(°>) (u)] ,T^aH-‘^)] iO) 



(3.75) 



which is a familiar expression for the polarizability in a CC framework [39]. We 
can rewrite Eq. (3.74) in the configuration basis [44] as 



= i:{0|(l+A(»')[p;,“>-{p™)]|h)A*q+a.)(h|pf-{Al”>)|0) + 






i:(0| (l + A'»l) [p;,“> - <p™)] |h)A-' (-a;) - <pi”'>|0> + 






(0| (l + A<»>) - E„) (i.,)tI'> (-C,) |0) 



(3.76) 
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where 



A (±w) = (h| ± uj) |h) , 



(3.77) 



where |h) indicates the manifold of singly and doubly excited determinants. = 

(0| (1 + A) and |f) indicates all determinants that are higher than twofold 

excited. In the limit that the |h) manifold is complete the quadratic term in Eq. 
(3.76) vanishes, or at the CCSD level the result is equivalent to Full Cl for two elec- 
trons. Alternatively, we can adopt the ’’propagator” viewpoint [21], which defines 



the polarizability as 



O^Aa (t^) = 




{'if - {/if )l'»f XiflA”’ - (Af 



/=0 k^O 






(3.78) 



(0) 



where and E)^\ A: = 0, 1, ... are eigenfunctions and eigenvalues of the unper- 
turbed zeroth-order Hamiltonian. The states form a complete biorthog- 

onal set (i.e. ') = 6ki). By using the CC parametrization for the ground 

and excited states (EOM-CC states), we can rewrite Eq. (3.78) as 



(t^) = E<0| (1 + A) I (pf - (Af >) |h)A-' ({-!)' O') {h|/if |0) (3.79) 

l-O 



which amounts to exactly the first two terms of Eq. (3.76). 

Equation (3.79) is also the same as Eq.(3.27) in the limiting static case, i.e. 
when a; = 0. We can therefore safely assume, that the time-dependent derivation 
leads to expressions which are equivalent to the time-independent formulas, while 
being much more laborous and complicated. 

However, extreme care has to be taken when going in oposite direction; from 
time-independent case to equations with frequency-dependence. Although general 
formulas remain the same, the presence of the frequency component causes the 
necessity to evaluate certain frequency dependent moments separately. The final 
property is not a sum of 24 identical moments as in case of second hyperpolariz- 
ability 7 [Eq.(3.29)j, but in general it becomes a sum of 24 different components. 
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each with its own structure of the frequency dependent resolvent. Of course, due 
to molecular symmetry, a great number of these components are identical and the 
general formulas for all the known electrooptic processes are well known [31, 45]. 

3.1.2 Role of the Renormalization Terms in the Sum-Over-States Expression 
Butcher’s formula for frequency-dependent polarizabilities in the notation of 
Bishop [32] is 






( ^cn • • • ^n) 




• • • 



5Z(^lA*al«l)(«l \l^l3 



Ci2 



fll 



a 



n 



■ ■ ■ {p-n \P'i \ 9} [(^ai ^cr){p^a2 ^cr 4" ^l) ‘ ’ ‘ ip^an ^n)] (3.80) 



where represents a summation over all forms obtained by simultaneously 

permuting pairs of dipole operators, Ji^ = fii — {g \fii \ g) and associated frequencies, 
uji. Also, LOfj = uj\ -\- UJ 2 -\- ■ ■ ■ + u)n. Writing this expression in terms of the resolvent 
operator, where |^>) = \g^ 02 . . . 



R{co,) = I#) {^\Eo-Ho + Lo, I#)-' (^1 = |<4>fc) ($fcl (uk - (3.81) 



^a,/3,...i = Pa,/9,...i PaR{-^a) PgR{oJl ~ UJ^)p^R{uJi + UJ 2 ~ ■ ■ ■ PiR{~UJn) Iq) ■ 

(3.82) 

In the EOM-CC theory, we have separate right {r} and left {1} biorthogonal 
eigenstates of if = e~^He^ that form a resolution of the identity, 1 = ]r)(l|r)“^(l[. 
For k = 0, (0| lo = (0| (1 -I- A), tq |0) = |0). Using the facts that |#) = |r) S = 

]h) C and ($] = S (1[ e~^ = C (h| where |h) = |0, fi, £ 2 , . . . f„) with ]0) the 
Fermi vacuum state and fi,f 2 , . . ., all single, double,. . . up to n-fold excitations, 
we can rewrite Fq. (3.81) as 



R{u;^) = e^]h)C(C“^(h F;o-iio + ‘^x h) ‘C-')C(h|e 



-1 



1-1 



-T 



(3.83) 



= e^|h)(h 



Eq — Hci E u) 



0 T U/j; 



h) (h| e ^ = e"* R(a;x)e 



T- 



-T 



7 
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with 



R-(^x) — (h £'o — Hq + ujx 



(3.84) 



Hence, a convenient, perfectly general form for the frequency-dependent hyperpo- 
larizabilities expressed in terms of unsymmetrical EOM-CC quantities is 



Y'n 



= Pa,/3,...i(0|(l •+• A)//a|h)R(-Wo.)(h|/i/3|h)R(a;i - o;^)(h|yu^|h)(3.85) 



R(cji -\- u>2 — a;(7-)(h|/ii-i|h)R(— a;i)(h|//i|0) , 



(3.86) 



0 



where /i* = e - (pj) , {ni) = (q (1 -h A)e 

The EOM-CC approach in its Cl approximation offers a convenient recursive treat- 
ment of second and higher-order properties. We simply consider the usual pertur- 
bation expressions that can be used to define such frequency-dependent properties. 
Hence, all the frequency-dependent polarizabilities can be written as Eq. (3.86). 
over the configurations |h). Since Eq. (3.86) is exact, it satisfies the time-reversal 
symmetry, and the exact form is independent of whether expectation values have 
the bra and ket states related as conjugates (see also [44]). 

Since we use EOM-CCSD, the first approximation is to restrict |h) to |0,fif2), 
which introduces a lack of extensivity, as in any truncated Cl approximation. 
But, it does not introduce any numerical time-reversal error even if it could for 
some approximations, as all computational expressions are built upon the average 
of frequency reversed expressions, as has been apparent from our first EOM-CC 
paper on polarizabilities [7], emphasized in refs [44] and [21], and is a customary 
expedient used, e.g. by Bishop [32]. Since including ]0) in ]h) causes singularities 
for static properties, we also consider its exclusion, or alternatively an expansion 
in the reduced resolvent i?o(<^x) = |fi, f 2 )(f|-S'o - Hq -h a;j;|f)~^(fi, £ 2 ], which leads 
to renormalization terms. We also need to address the issue of renormalization 
terms. Erom consideration of the inverse, Eq. (3.84) can be readily expanded to 
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give the relation 

R{<^x) = ^ |0) (0| (1 + ~ |0) (0| (1 + ^)-Ro(<^x) + Ro{<^x) 7 (3.87) 

where Rq{ux) = |f)Ro(a;a;) (f|. Substituting Eq. (3.87) into Eq. (3.82) we have 

^a/?(~^o-7 ^l) = Pa/3 (0 1(1 + A)flaRoi‘^l)Ji0\ 0) (3.88) 

since the other two terms must vanish by virtue of (0 |(1 + A)/1 q| 0) =0. Then for 
the second- and third-order susceptibilities, 

^^l7 ^2 ) = Pq/ 37[(0 |(l -f- A)//q7?o(‘^1 + ^2)M/9-Ro (^ 2)^7 1 0) 

— (0 1(1 -I- A)licRo{uJi + ta2)/l/3| 0) (0 1(1 -I- A)/?o(^ 2 )/l 7 l 0)] 

(3.89) 

Pq/37(S[(0 1(1 + A)/J,aRo(Ui -|- 032 -|- U 3 ) Hi3Ro{u>2 + 0;3 ) //.y i?o (^3 ) 1 0) 

— (0 1(1 -I- A)JlaRo{uJi + LO 2 + L 03 )]lpRo{uJi)liy\ 0) (0 |(1 -h A)7 ?o(^^3)M(5| 0) 

— (0 1(1 -I- A)JiaRo{uJi LO 2 + u>z)]i.0\ 0) (0 1(1 -h A)RQ{uJi)Jl^R(j{LO^)lii\ 0) 

— (0 1(1 -|- A)/iccRo(<^i + a ;2 + o; 3 )/x/j| 0) (0 |(1 -I- A)/l^i?o(^i)-Ro(‘^3)/l(5| 0)] . 

(3.90) 

which is a simple derivation of polarizability expressions without secular diver- 
gences [32]. The renormalization terms are exactly what one would expect from 
elementary considerations based upon inhomogeneous equations. For ca 7 ^ 0, an 
equally good computational procedure is to simply use Eq. (3.86) with jh) = 
| 07 fi 7 f 2 ) which avoids renormalization terms. 

Once the configuration space is restricted the neglected higher excitations have 
two roles; to add additional correlation corrections, and, to cancel the unlinked 
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Table 3.1: Role of the renormalization terms at equilibrium geometries. 



system / basis 


mr 


mr mr 


7(ll)“ 


7(11)" 


7(11)^ 


7(11)'^ 


HF/POLl 


-5.44 


-6.17 -5.76 - 6.08 


370 


406 


391 


403 


HF/POL1++ 


-7.10 


-8.03 -7.67 - 7.89 


466 


511 


499 


504 


N2/POLI 


- 


- - - 


950 


996 


978 


991 


N2/POLI++ 


- 


- - - 


990 


1027 


1009 


1024 


CO/POLl 


23.0 


24.5 23.8 24.1 


1480 


1519 


1494 


1515 


H2O/POLI 


-12.1 


-14.2 -13.4 - 13.9 


1210 


1308 


1269 


1278 


H2O/POLI++ 


-14.9 


-16.4 -15.4 - 16.0 


1710 


1792 


1744 


1768 


C4H6/ 

6-31G+(P,D) 




— — — 


18480 


20230 


19610 


20030 



“ finite field derivative of the orbital relaxed dipole moment 
** EOM-CCSD analytical value without renormalization terms; EOM(l) 

^ EOM-CCSD analytical value with renormalization terms; EOM(2) 

^ finite field derivative of the EOM-CCSD quadratic polarizability (exact CCSD 
unrelaxed value) 



diagrams (extensivity error) left in any truncated approximation. The second can 
be handled by including certain quadratic, cubic and quartic terms in 
and respectively, which arise from coupled-cluster energy derivatives that 

correspond to certain higher-excitation effects [44, 46]. To regain the formally 
attractive perturbation theory form, which is only possible within configuration 
space, plus its convenient recursive generalization to any order, we prefer not to 
include the computationally difficult cubic and quartic terms, so we depend upon 
calibration to the available exact results in a given basis to offer confidence in 
the values. Eor static hyperpolarizabilities, such calibration is accomplished by 
comparing with finite field derivatives of the energy, dipole, or polarizability in the 
presence of a field. The procedure can also be used for the Kerr effect [47]. For 
other frequency-dependent results, we also compare to the known results [48]. All 
of these were employed to help justify our approximate results for the HE molecule. 
In Table 3.1, we show such an analysis. 
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Table 3.2: Percentage error of the EOM-CCSD static hyperpolarizabilities for 
(LiH)„ 

1.59725 a.u., 3-21G basis set. All values in percent, unless noted otherwise. 



number of LiH 


/j( )“ 


/?( )" 


/ 5 ( 


1 


1.2 


0.6 


-501 a.u. 


2 


2.0 


0.8 




3 


3.1 


1.1 




4 


4.0 


1.5 




5 


5.4 


1.9 






7( )“ 


7( )‘ 


7( 


1 


0.4 


0.1 


78 580 a.u. 


2 


0.5 


0.1 




3 


0.6 


0.1 




4 


0.7 


0.1 




5 


0.7 


0.1 





“ EOM-CCSD Cl-like without renormalization terms; EOM(l) 

* EOM-CCSD Cl-like with renormalization terms; EOM(2) 
finite field derivative of the dipole moment. 

Neglecting the renormalization terms gives values which are actually closer 
to the exact values than if they are retained. Furthermore, the effect of orbital 
relaxation shown in the first column for j3 and 7, is far more numerically important 
than the differences in either of the EOM approximations, which, here, always 
bracket the exact value. It is evident, that either EOM approximation gives results 
within the estimated error bars of the calculations. For HF/POLU— I- /?(||) the 
orbital relaxation effect is -10%, while EOM(l)’s error is -1-1.8%, EOM(2) -2.8%, 
and the triples estimated effect is 5%. 

We have also considered the problem of size extensivity of the resulting ex- 
pressions (Table 3.2). Given a chain of non-interacting LiH molecules, in a 3-21G 
basis (all electrons), the percentage error for the first hyperpolarizability /3(||) at 
LO = 0.0, calculated using the unrenormalized EOM(l) method is 1.2% for one 
molecule and 5.4% for 5 molecules. The same errors computed with the EOM(2) 
approximation are 0.6% and 1.9% respectively. Percentage errors for the second 
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hyperpolarizability 7(H) at u> — 0.0 are 0.4% and 0.7% for 1 and 5 LiH molecules 
in EOM(l) and 0.1% for both 1 and 5 LiH molecules in EOM(2). 

Apparently, the renormalized EOM(2), which is generally in slightly poorer 
agreement with reference results, has a somewhat smaller extensivity error (al- 
though that does not mean it is closer to the complete Cl limit). Eor trans- 
butadiene in the POLl-l— I- basis, the differences between EOM(l) and EOM(2) 
for 7(H) are 0.72 at uj = 0.0 a.u., 0.84 at to = 0.0430 a.u. and 1.18 at a; = 0.0656 
a.u. at the experimental geometry, compared to the estimated error bars of ±2 
a.u. in the calculations. 

Since there is no orbital relaxation in the full Cl limit, the remaining relaxation 
effect is partly indicative of how far from that limit the calculations remain. 

3.2 Brillouin-Wigner Perturbation Theory Based Methods 

In modern quantum chemistry of nondegenerate systems in which the ground 
state is well described by a single-determinant wavefunction, the Rayleigh-Schrodinger 
Perturbation Theory based methods (MBPT and CC) are very efficient for the de- 
scription of properties. The most attractive feature of these methods is their size 
extensivity (proper scaling with the number of particles). Since size extensivity is 
a very important property of a method, it is expected that all the newly developed 
methods will satisfy this property. Possibly this is the reason why less attention 
has been paid to Brillouin-Wigner Perturbation Theory [49, 50, 51, 52]. Redefining 
A in the perturbation expansion 3.4 gives us the opportunity to shape the form 
of the wavefunction and energy expansions. Choosing A to be the total energy of 
the perturbed system E, we obtain the Brillouin- Wigner expansion, giving rise to 
BWPT. Putting X = E cancels —E -|- A in the resolvent expression 3.6 leads to 

00 

m=0 



(3.91) 
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for the wavefunction and 



A£; = i^o\v {ME)vf 

m=0 




(3.92) 



for the energy expression. 

The expansion in zero-order functions gives the following terms: 

(^Hq — (first — order equation), 

[Ho - (^(1) _ |$(i)) + ^(2)|$(o)^ (second - order), 



(3.93) 

(3.94) 



[Ho - = (^(1) - V) |$(2)) + + ^^^’|$(1)) 

(third — order), (3.95) 

[Ho - E^^^) - V) + E^^^l^il)) + 

£;(2)|$(2)) (fourth — order) . (3.96) 



After applying the left coupled-cluster ground state eigenvector and express the 
perturbed wavefunctions in terms of the zero order CC wavefunction we obtain 



($o| (Ho - E^^^) |<J>(i)) = ($o| {e^'^ - V) |$o) , 



(3.97) 



(«4>o| (^0 - E^^'>) |$(2)) = ($q| |<j,(D) + ^( 2 )|,j>(o)^ ^ 

($o| {e^^^ - V) R{E)oV\^o) + , (3.98) 

($o| (^0 - E^°^) = ($o| {e^^^ - V) |<I>(2)) + f;(3) + £;(2 )($o|${i)) 

($o| {e^^^ - V) R{E)oVR{E)oV\^o) + E^^^ + E^^\^o\R{E)oV\%) , 

(3.99) 



($o| (^0 - E^^^) = ($o| (e^^^ - V) + E^^^ + 
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($o| - V) RiE)oVR{E)oVRiE)oV\%) + E^^^ + E^^\^o\R{E)oV\%) + 

E^^^^o\R{E)oVR{E)oV\^o) • (3.100) 

In above equations R{E)o is the BW resolvent operator 

R{E)o = 

ij 

= ^ (3.101) 

ij 

defined for a non-diagonal H Hamiltonian. 

The perturbed wavefunction |$(n)) can be written in terms of the unperturbed 
ground state wavefunction with the excitation operators being actually perturbed: 

|${3))^^ — J"(3) _j_ j'(2)j’(l) _|_ 7"(0j'(l)j'(l)|<^p\ 

= J'(4) _j_ 2^(3) j>(l) _|_ j>(2)ji(2) _|_ J'(l) j'(l)j’(l)2^(l)|^^\ 

(3.102) 



Applying the above results to the zero order expansions in Eqns. (3.93) to (3.96) 
and noting that (d>o| [Hq — T(”^|$o) = 0 we get 



E(^) = (d>o|I>|$o) , 



E(^^ = (d>o| (e(^^ - V) R(F)oVl^o) + ($o| (^0 - E(°^) , 



(3.103) 

(3.104) 



= ($o|(^‘'^-t>)i?(^)oT/?(£;)oI>|$o) 



-h 






+ ($o| (^0 - E(^^) T(2)T(i)|d>o) , 



(3.105) 



E(^^ = ($o| - T) R(E)oVR(B)oVR(E)oVI^o) ' 



1 
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+ + E^^^^o\R{E)oVR{E)oV\^o) 



+ 



+ 



+ 



+ 



($ol (^0 - E^°^) 

($ o | (^0 - 
($o| {Ho~ 

(<4>o| . 



(3.106) 



In the same way as with the RSPT equations (3.27) to (3.29) we can derive 
expressions for static BW-EOM-CC polarizability and hyperpolarizabilities. Inter- 
estingly, the higher order terms coming from expanding are exactly the same 
as in RSPT and they are completely independent on the choice of the resolvent 
and therefore time and frequency independent. Both this reason, the fact that 
the numerical significance of such expressions has been proven to be marginal [53] 
and difficulty and cost of computational implementation caused us to neglect these 
terms. 



a(0; 0) = 2(0|(1 + A)RRo(a)P|0) , (3.107) 



/?(0;0,0) = 6(0|(l + A)I/Ro(/?)Ei?o(^)P|0) + 
• 3a(0;0)(0|(l + A)Ro(/5)P|0) , 



(3.108) 



7(0; 0,0,0) = 24(0|(1 + A)RRo(t)^^o(7)^^o(7)^|0) + 

3/3(0; 0, 0)(0|(1 + A)Ro(7)^|0) + 6«(0; 0)(0|(1 + A)Ro(7)^^o(7)V^|0) • 

(3.109) 

The equations (3.107) to (3.109) differ from similar RSPT equations by having 
fewer renormalization terms (had the Hamiltonian be symmetric and hermitian, no 
renormalization terms would be present in this approximation), by replacing the IT 
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operator by V (no extra unlinked terms coming from — (0|V|0)) and by including 
the exact, unknown energy on both sides of the equation. The only plausible 
method of solving this type of equations is by iterating; the starting value of E (in 
most cases it is the second-order RSPT perturbed energy which is relatively easy 
to determine) the right hand side of Eqns. (3.107) to (3.109) is solved. The newly 
determined value of E is then inserted in the right hand side of these equations and 



accuracy. This 



the process is repeated until the value of E converges to desired 
feature, along with the slow convergence of the series and lack of size extensivity 
caused the BWPT method to be unpopular among quantum chemists, although 
it has been used extensively in nuclear physics [33]. In quantum chemistry the 
BWPT was first applied by Lowdin [54, 55, 56]. Related to this problem was 
the works of Brandas and Bartlett [57, 58] and more recently by Hubac et al 
[59, 60, 61]. Thanks to the pioneering work of Hubac [61] we know that in the 
single-reference case, the BW theory based on the coupled-cluster wavefunction 
is fully equivalent to the standard CC theory and although the BWCC does not 
employ the Baker-Campbell-Hausdorff formula it is a fully size-extensive method 
since the iterative procedure used in solving the BWCC equations cancels out all 
disconnected diagrams (see Appendix 1). Because of the choice of our Hamiltonian, 
the right and left eigenvectors are not each other conjugates-they are biorthogonal. 
Thus, our left eigenvector expressed in terms of the ground state (Fermi vacuum) is 
($o| = (0|(1 + ^) where A is a linear, Cl-like deexcitation operator, complementary 
to the T excitation operator [42]. The right eigenstate is defined as l^o) = |0). 

Since the equations for static processes are readily extendible to evaluate the 
dynamic, frequency dependent quantities (as shown in Sect. 3. 1.1) and starting with 
the Fermi vacuum ground state, the final equations for (hyper)polarizabilities in 
the BWPT based EOM-CC are as follows: 



a^j{-u-u) = (0|(1 + A) -H,)i?(a;,a)H,-|0) 
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+ (0|(l + A)(£;(i)-1/^)i?(-a;,tt)V;|0) (3.110) 

(3ijk{-uJa-,u;uUJ2) = Pihuj) ((0|(1 + A) - t>,) R{co„ P)VjR{u2, mk\0)) 

+ P(i, oj) ^o;jj(t(;cr)(0|(l + A)i?(ci;2, /?)^fc|0) j (3.111) 



7u'A:/(~^t7; ^25 ^3) — 

Pii^cu) ((0|(1 + A) (e(') - k) R{cJa^l)VjR{uj,+U2n)V^^^ 

+P(z, u) C(;2)(0|(l + A)/?(cl;i, 7)t^|0)^ 

+P(z, u) ^a2j(~^3)(0|(l + !s)R{lo2 + ^1, 7)t4P(c^^i, 7)^10)) (3.112) 

where u>a is a sum of all the incoming frequencies, u;i,...,c<;„ and is replaced 
by o;, /?, 7 in the resolvent 3.101. 

As mentioned earlier, the Brillouin-Wigner resolvent R{E) is no longer secular 
at the excitation energy frequency: while the RSPT resolvent goes to 0 when 
the incoming radiation has the frequency corresponding to an internal electronic 
excitation {lo = E — Hq), the BWPT resolvent remains finite due to a damping 
from the extra energy factor Ei (3.101). This enables us to study the nonlinear 
processes near and at the resonant frequencies in a fully consistent way. 

The current progress in theoretical non-linear optics involves mostly non-resonant 
processes: there exist several models for calculation of hyperpolarizabilities, both 
optical, magnetic and magneto-optic. All of them however contain singularities 
at the poles of electronic excitations, making it impossible to determine the ac- 
tual magnitude of a process at the resonance frequency. The divergence due to 
resonance is clearly unphysical, and ultimately has its origin in the neglect of 
higher-order terms and the quantum field interaction. When these are taken into 
account, large changes in the electron density caused by the resonance introduce 
level shifts. The general and widely accepted description of the non-linear effects is 
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based on the anharmonic oscillator model (Sum-Over-States expression) contains 
a damping factor in the denominator. This damping factor plays role of the level 
shifter here, but the actual choice of this quantity is based on empirical observa- 
tions like natural lifetime in absence of collisions etc. This factor can be eventually 
obtained analytically by summation of the higher order terms [62]. 

Another way of a level-shifting technique is the Brillouin-Wigner perturbation 
expansion in which one has an exact energy on both sides of the SOS equation and 
the price paid for the flexibility of the method is that these equations have to be 
solved iteratively. 

For more complicated properties with many different incoming and outgoing 
frequencies every denominator should in principle be the Brillouin-Wigner type 
level-shifted denominator. Since the energy shifts are also direction dependent 
(they are just the processes being computed), the convergence in such a case is 
very slow and the cost of such calculations is prohibitive. 

A solution to this dillema is quite simple: because at any given frequency there 
can be only one resonant denominator, we can apply the BW level shifting to it 
and leave all the other denominators out of the iteration procedure, assigning them 
a standard shift (static value of the property in a given direction). 

The numerical effect of such an approximation is expected to be very small as 
the leading, resonant term is much more significant numerically (typically 2-3 or- 
ders of magnitude) than the nonresonant terms. The example of these approxima- 
tions is presented in Table 3.3. Here we compare the resonant linear polarizabilities 
of the HF molecule in the DZP polarized basis set. 

The first polarizability is chosen for this example because the expression con- 
sists of only two terms and the explanation is rather straightforward. Obviously, 
the RS-EOM-CC method fails at the reonance due to the structure of its denomi- 
nator (first column of the table). However, applying Eq.3.110 removes zero value 
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Table 3.3: Resonant linear polarizabilities of the HF molecule. 



process/method 


RS-EOM-CC 


BW-EOM-CC“ 


BW-EOM-CC* 


0!zz{+^) 


inf 


2699.3 


2695.8 




1.1348 


1.1348 


1.1346 



“ Full form; from Eq. 3.110 
** Approximated form; from Eq. 3.113 

from the SOS denominator and delivers a reasonable value of the first polarizability 
(second column, BW-EO-MCC full). Of the two terms in Eq. 3.110 only one is at 
the resonance for either positive or negative value of the external electric field. Ap- 
plying the full BW-EOM-CC expression is quite expensive computationally, as it 
requires 14 iterations to fully converge. However, if one substitutes the exact value 
of a in the non-resonant term by a constant (static polarizability value) and iter- 
ate the equation changing only the resonant term polarizability, the convergence 
is much faster. The resulting equation 

= (0|(1 + A) (iF(i) - t>) R(-fa;,a(0))l>|0) 

+(0|(1 + A) (^(1) - K) R(-tj,a(a;))lA|0) (3.113) 

converges in 6 cycles and delivers almost the same value of the first polarizability. 

For higher-order processes we may expect much higher computational cost in 
converging the BW-EOM-CC equations if all the frequency dependent moments 
are of the BW type. 



CHAPTER 4 

APPLICATIONS TO MOLECULAR SYSTEMS 



4.1 Eirst Polarizability and Long-Range Dispersion Coefficients 
The working equation for EOM-CCSD polarizability has been already derived 
(Eq. 3.79) and it reads 

= 1](0| (1 + A) I (/ij’ - (/ij')) |h)A“^ ((-!)' w) (h|/ij^|0) , (4.1) 

1=0 

where A is defined as 

A (±o;) = (h| ± to) |h) . (4.2) 

This provides a computationally and conceptually convenient alternative approach 
to polarizabilities. The EOM-CC has the advantage of providing the and 

E'fc \ A: = 0, 1, ... in Eq. (3.78). From the above development it is clear that this 
is not the full derivative for a truncated CC method, but it offers an equally valid 
’’propagator” viewpoint that has often been used in SOPPA and CCPPA property 
calculations [63]. The analogy with an expectation value and the first-derivative of 
the energy relative to a field is apparent. If the (generalized) Hellmann-Feynman 
[64] theorem is satisfied, the two forms are identical. If not, different results are ob- 
tained. Depending upon your viewpoint, either can be the ’’rigorous” definition of 
the property. The same condition applies for second- and higher-order properties. 
If the higher analogue of the Hellmann-Feynman theorem is satisfied, it means that 
the results from ordinary perturbation theory and derivative theory are identical. 
Since the derivative always introduces the appropriate non Hellmann-Feynman 
terms, it should be somewhat superior numerically. The straight ’’propagator” 
approach using the EOM-CC states is referred to as the Cl-like approximation 
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[21], as the excited states in EOM-CC are obtained from a Cl-like diagonaliza- 
tion procedure of a transformed Hamiltonian based upon a CC ground state. The 
EOM-CC approximation for excitation energies is not entirely linked [65], like Cl. 
Second order properties in EOM-CC are fully linked in the quadratic approxima- 
tion but not in the Cl-like approach. This scarcely affects a single molecule, but if 
we replicate the molecule many times we would find a numerical problem with the 
Cl-like approximation in the limit [66]. We can correct this, however, by removing 
the unlinked terms that remain in the Cl-like approximation. Eormally this means 
we take from the quadratic term that needed to correct the unlinked diagrams in 
the lead term, and then in the interest of efficiency, neglect the remaining linked 
quadratic part. This defines a linear linked approximation [67, 44] 

= i:(0| (1 + A') |p!,Vjh>A-l((-l)'a;))(h|/l<»»|0) . (4.3) 

1=0 



Here 



_( 0 ) -T T 

1^ A, open ^ 



(0|e V^e^lO) , 




while A' is an explicitly connected operator that is defined from a slightly modified 
A equation obtained by equating the usually small matrix elements 

{0\e-'^He'^O = F„ + ('1-5) 

responsible for disconnected contributions in A, rigorously to zero. The linear 
approximation scales properly for any number of molecules, but is no longer ’exact’ 
for two electron systems. However, it retains the convenient ’’propagator” form, 
shared by the exact result. In practice we replace A' by A since this changes the 

t 

results only very slightly, but it allows us to calculate both the Cl-like and the linear 
approximations in a single calculation. It follows that in the present calculations. 
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if a molecule does not have a permanent dipole moment (due to symmetry) the 
linear and Cl-like approximations yield identical results. In the next section we 
will also consider the fully linked linear approximation, and we will refer to this as 
linear’ (in accordance with A'). Equations (3. 76), (3. 79) and (4.3) (with A' = A) 
summarizes the quadratic, Cl-like and linear approximations to polarizabilities. In 
the following we will consider all three approximations numerically. 

The difference between the Cl-like and linear approximation to the polarizabil- 
ity can succintly be put as 

/ \CI —like / \ Linear 

(^^) = CtAa 

+ ((0|/i^|0) - (0| (1 + A) A^IO)) (0| (1 + A) • Ta (u;) |0) 

+ ((0|/2„|0) - (0| (1 + A) flam (0| (1 + A) • Ta (-ui) |0) . 

(4.6) 

Therefore, in general, if the dipole moment of the reference determinant and the 
correlated CC dipole moment are similar, there is very little difference between the 
linear and Cl-like approximations. For very large systems there is a problem with 
the Cl-like approximation, however, because the term involving the perturbed 
amplitudes can grow indefinitely. For larger systems the linear approximation, 
therefore, offers distinct advantages over the Cl-like approximation [67]. Most 
of our calculations are limited to small to medium sized systems and the differ- 
ence between the linear and Cl-like approximation due to improper scaling of the 
Cl-like approximation is usually negligible. The other limiting case occurs when 
there is a large difference between the correlated and Hartree-Fock dipole mo- 
ments. As an example we consider the HF molecule at three internuclear distances 
Ro=T7328, 3.0 and 10.0 a.u. At the RHF level the system incorrectly separates 

f 

into -I- F~ fragments, and, correspondingly, carries a large dipole moment. 
The CCSD calculation brings this back to the proper ground state, separating into 
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neutral fragments, H + F, which in the limit does not have a dipole moment. This 
is ilustrated in table 4.1, where we list the dipole moments for the three distances 
considered. In table 4.1 we also include the CCSD amplitudes. These amplitudes 
are very large at larger separation indicating again that the Hartree-Fock refer- 
ence state is not a good description of the true ground state. If we subsequently 
calculate the static polarizabilities we obtain the results presented in table 4.2. It 
is obvious that the Cl-like and quadratic results agree quite well at all distances, 
but the linear version does not follow this pattern, and is in fact erroneous. Dia- 
grammatic analysis shows that in the linear approximation, so-called EPV terms 
are neglected, and the present pathological example highlights these contributions. 
The situation is similar to the correlation energy. For two-electron systems CCSD 
and CISD are exact, and both can be said to include EPV terms. On the other 
hand, linearized CCSD does not include such contributions, but unlike CISD, it 
is rigorously size extensive. However, it is not exact for two-electron systems. In 
fact LCC(S)D usually overshoots the correlation energy. The situation for second 
order properties is a little different, because the extra (EPV-related) term in the 
Cl-like approximation is usually very small. We have also included the fully linked 
linear’ results. Compared to linear these results deviate further from the more 
correct quadratic results for larger separations. This can be understood because 
the disconnected contribution to A is quite appreciable (it follows the pattern for 
t\) , and in the linear’ approximation we neglect such contributions from A, even 
though the disconnected contributions in A mainly give rise to connected contri- 
butions to the polarizability. In this context the extended coupled-cluster method 
(for analysis see ref. [64]), which treats both the right and left hand ground state 
in exponential fashion is most satisfactory. 
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Table 4.1: Dipole moments of the HF molecule. 

Rh-f (a.u.) /i (SCF) n (CCSD) ti coefficients 
1.7328 (Re) 0.75628 0.70655 0.021 

3.90 1.78678 0.58405 0.182 

10.0 4.03577 -0.07291 0.396 

Table 4.2: a^z component of the static polarizability for the HF molecule. 

Rh-f (a.u.) Cl-like Linear Linear’ Quadratic 

1.7328 (Re) 6.568 6.570 6.590 6.498 

3.90 28.49 31.38 31.63 28.48 

10.0 8.461 9.922 11.21 8.410 



We conclude that in the single molecule calculations that are of practical in- 
terest to us, the Cl-like approximation does not suffer significantly from the size- 
extensivity error. Due to the inclusion of EPV terms we actually expect it to be 
a little more accurate than the extensive linear approximation for some cases of 
interest. In all of the examples considered in this thesis (except the above patholog- 
ical example) the Cl-like and linear approximation yield identical results up to the 
figures quoted in this thesis. This is true exactly if molecules do not have a dipole 
moment due to symmetry. Let us emphasize finally that the EOM quadratic model 
as the derivative is formally most satisfactory, but it looses the convenient propa- 
gator form and can become expensive computationally. This is particularly true if 
we apply partitioning techniques to the EOM scheme [68, 69, 70]. We know this 
is a valid approach to polarizabilities as well as NMR spin-spin coupling constants 
[70], and in the partitioned model we definitely do not want to use the quadratic 
model, since it would forfeit all time-savings obtained by the partitioning. 

In the following we wish to make a survey of several molecules in a consistent 
basis, rather tham the ulimate converged result for one example. Hence, we use ex- 
perimental geometries and Sadlej’s polarizability basis, POLl [71, 72] which consist 



52 



of a (14s,10p,4d)/[7s,5p,2d] contraction for S, a (10s,6p,4d)/[5s,3p,2d] contraction 
for C, N and 0 and a (6s,4p)/[3s,2p] contraction for hydrogen. For trans-butadiene, 
to be consistent with prior work [5, 73], we have used the standard 6-31G basis set 
augmented by two diffuse p and d shells with the same exponent Cp = Cd - 0.05 
which has been previously shown to combine small size with good accuracy for 
this molecule. Cartesian Gaussian basis functions have been used in all the cal- 
culations. All the results reported were obtained by using the ACES II program 
system [74]. 

The dipole polarizability is a diagonal second-rank tensor. For any linear 
molecule two of these three components are identical due to symmetry. The two 
unique components are commonly referred to as the perpendicular and par- 
allel a\\ components with respect to the principal rotational axis. The average 
polarizability < cr > and the polarizability anisotropy Acr are the quantities most 
commonly determined experimentally [75]. They are defined as 



< Q: > 

Aa 



a (xx) 4- a (yy) 4- a {zz) 



3 

[(« (xx) - a {yy)f 4- {a (xx) - a {zz)f + (a (yy) - a {zz)fYl‘^ 




which reduces to the simple expressions for molecules possessing a 3-fold or higher 
rotation axis (e.g. a linear molecule). 



2q;x 4- ftii 





All the molecules considered in this study are linear, and hence we report the 
average polarizability and the polarizability anisotropy (henceforth referred to as 
< a > and Ao respectively. The Ce dispersion coefficients are calculated from 
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polarizabilities evaluated at selected imaginary frequencies ia using 

<^6 = 1 ^ 1 / {t) dt 



and Gauss-Legendre quadrature where the function ^(t) is defined as 







(4.10) 



by means of the substitution 




(4.11) 



The parameter uo has been chosen to be 0.1 a.u. as in Amos et al. [76]. 

Nitrogen, Carbon Monoxide and Acetylene. 

In tables 4. 3, 4. 4 and 4.5 we present the calculated dipole polarizabilities of 
the N 2 , CO and C 2 H 2 molecules at several different frequencies along with the 
corresponding experimental results. The ‘CCSD‘ results are obtained by the finite- 
field method as the difference of the analytically computed dipole moment [77] 
with orbital relaxation, and are limited to the static values. The TDHF values are 
from the analytical derivative program of Sekino and Bartlett [78]. For N 2 and 
CO we have also included second-order polarization propagator approximation 
(SOPPA) results [79], while for N 2 , there are MC-TDHF results. Both the latter 
are in different basis sets from the CC results but agree exceptionally well with the 
EOM-CC results. As we can see from Tables 4.3 to 4.5, the electron correlation 
effects in general are fairly minor for these molecules, the TDHF and correlated 
dispersion for the isotropic component of the polarizability agrees fairly well for 
all methods. However, the most notable exception is the dispersion for Aa of the 
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CO molecule which is negative (but very small) at the TDHF level, and positive 
(but small) at the correlated level. This is likely to be related to the well-known 
sensitivity of the (small) dipole moment of CO which changes sign upon inclusion 
of correlation. In the case of acetylene the agreement is less satisfactory, even for 
the isotropic parts, while also the discrepancy between TDHF, SOPPA and the 
EOM models is more pronounced. 

Orbital relaxation effects are very minor as follows from the comparison be- 
tween finite-field-CCSD and the EOM-CC quadratic model for all three molecules. 
It is well known that CCSD shows remarkable insensitivity to orbital choice [80, 43], 
but the EOM-CCSD (or CCSD linear response) will not necessarily be as insensi- 
tive. Here, however, the ground state molecular orbital relaxation effects are very 
small. Finally the static, absolute, values for the Cl-like and quadratic EOM-CC 
models agree quite nicely. Inclusion of correlation improves the agreement with 
experiment, but TDHF results are already pretty good. 

Chlorine 

The calculated dipole polarizabilities of the CD molecule for several different 
frequencies, and the corresponding experimental values are given in Table 4.6. 

Here there is a substantial difference between the correlated dispersion and that 
given by TDHF. The EOM Cl-like and quadratic model dispersion curves agree 
very well as before and there is only a minor effect from orbital relaxation. How- 
ever, correlation effects for CI 2 are more significant. For example, the correlation 
contribution to the static polarizability is 3.10 a.u. (10% of the total correlated 
(finite field CCSD) result) for \ai and 1.57 a.u. (0.94% of the total correlated re- 
sult) for Ao;. The EOM-CCSD Cl-like results overestimate the electron correlation 
effects by almost 1 a.u. compared to experiment, whereas the quadratic results 
lower the deviation to about 0.5 a.u. The role of correlation (and/or basis set) 
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^ G. R. Alms et al. [82]. 

® Bridge and Buckingham [83]. 
^ Luo and Jorgensen [84]. 
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Table 4.4: Dynamic polarizabilities of the CO molecule“. 



uj (a.u.) 


TDHF'’ 


SOPPA'^ 


CCSD'’ 


EOM-CCSD" 

Cl-like Quadratic 


exp. 


.0000 


12.23 


12.45 


12.95 


(c(0)) 

13.28 


13.07 


00 

o 

CO 


.0720 


0.22 


0.24 




{a{uj)) - (a{0)) 

0.27 


0.25 


0.27^^ 


.0886 


0.34 


0.37 




0.41 


0.39 


0.41'' 


.0934 


0.38 


0.41 




0.46 


0.43 


0.47^^ 


.0995 


0.44 


0.47 




0.52 


0.50 


O.OO'' 


.0000 


3.37 


4.45 


3.97 


Aa(0) 

4.39 


4.17 


3.59 


.0720 


-0.01 


0.07 




Aa{uj) — Aa(0) 
0.03 


0.02 


3.59^ 


.0886 


-0.01 


0.10 




0.03 


0.03 


- 


.0934 


-0.02 


0.11 




0.03 


0.04 


- 


.0995 


-0.02 


0.12 




0.04 


0.04 


- 



“ At the experimental equilibrium bond length R(CO) = 2.132 a.u. 
^ Present work 

Oddershede and Svendsen [79]. 

International critical tables [81]. 

® Extrapolated from International critical tables [81]. 

^ Bridge and Buckingham [83]. 



Table 4.5: Dynamic polarizabilities of the C 2 H 2 molecule“. 



Lu (a.u.) 


TDHF'' 


GGSD'’ 


EOM-CCSD‘' 

Cl-like Quadratic 


exp. 


.0000 


23.11 


22.17 


(a(0)) 

22.60 


22.25 


22.67'' 


.0720 


0.65 




{a{uj)) - {a{0)) 
0.58 


0.57 


0.80'= 


.0886 


1.00 




0.90 


0.89 


l.OR 


.0934 


1.12 




1.01 


0.99 


1.28^ 


.0995 


1.28 




1.15 


1.13 


1.4R 


.0000 


12.52 


12.30 


Aa(0) 

12.97 


12.50 


11.84'' 


.0720 


0.34 




Aa{uj) — Aq;(0) 
0.31 


0.30 


0.54^ 


.0886 


0.52 




0.48 


0.46 


0.74^ 


.0934 


0.58 




0.53 


0.51 


0.88" 


.0995 


0.66 




0.61 


0.58 


1.08" 



“ The equilibrium experimental bond length r(CH) = 2.004 a.u. and r(CC) 
2.274 a.u. 

** Present work. 

G. R. Almset al. [82]. 

Extrapolated from G. R. Alms et al. [82]. 
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Table 4.6: Dynamic polarizabilities of the CD molecule". 











EOM-CCSD" 




Lo (a.u.) 


TDHF'' 


SOPPA" 


CCSD'' 


Cl-like 


Quadratic 


exp. 


.0000 


27.38 


24.42 


30.48 


(a(0)) 

31.32 


30.86 


30.35" 


.0720 


0.84 


0.47 




(a(a;)) - (a(0)) 
0.66 


0.66 


0.72'' 


.0886 


1.31 


0.74 




1.04 


1.03 


1.11'' 


.0934 


1.48 


0.84 




1.18 


1.15 


1.23'' 


.0995 


1.70 


0.98 




1.37 


1.32 


1.41" 


.0000 


15.22 


24.44 


16.79 


Aa(0) 

17.46 


16.94 


17.56-^ 


.0720 


0.95 


0.86 




Aa{uj) — Aa(0) 
0.58 


0.56 




.0886 


1.50 


1.29 




0.91 


0.85 


- 


.0934 


1.69 


1.43 




1.04 


0.95 


- 


.0995 


1.96 


1.60 




1.20 


1.07 


- 



“ At the equilibrium experimental bond length r(ClCl) = 3.755 a.u. 
** Present work. 

" Oddershede and Svendsen [79]. 

" International critical tables [81]. 

" Extrapolated from International critical tables [81]. 

^ Bridge and Buckingham [83]. 



effects is further demonstrated by the rather poor agreement of SOPPA results 
with the experimental values. 

Carbon Dioxide, Carbon Oxysulfide and Carbon Disulfide. 

In Tables 4.7 to 4.9 we present calculated < cr > and Aa values for the CO 2 , 
OCS and CS 2 molecules for several different frequencies. The results obtained by 
the TDHF method and experiment are also given for comparison. 



It is evident from tables 4. 7, 4. 8 and 4.9 that orbital relaxation effects are very 
important for these molecules, and this importance increases with the number of 
sulfurs. The orbital relaxation effect (EOM quadratic - CCSD) for the isotropic 
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Table 4.7: Dynamic polarizabilities of the CO 2 molecule". 









EOM-CCSD" 




u) (a.u.) 


TDHF'’ 


CCSD'’ 


Cl-like 


Quadratic 


exp. 


.0000 


15.84 


17.44 


(a(0)) 

18.47 


18.00 


17.52" 


.0720 


0.20 




(tt(cj)) - (o;(0)) 

0.29 


0.27 


0.27“^ 


.0886 


0.30 




0.44 


0.42 


0.34" 


.0934 


0.33 




0.49 


0.46 


0.47" 


.0995 


0.38 




0.56 


0.53 


0.50" 


.0000 


12.08 


14.43 


Aa(0) 

16.32 


15.40 


13.73" 


.0720 


0.26 




Aa(a') — Aa(0) 
0.41 


0.40 


0.54<^ 


.0886 


0.39 




0.63 


0.60 


0.68" 


.0934 


0.44 




0.71 


0.68 


0.88" 


.0995 


0.50 




0.80 


0.77 


0.95" 



“ The experimental equlibrium bond length r(CO) = 2.195 a.u. 
* Present work. 

" G. R. Alms et al. [82]. 

Extrapolated from G. R. Alms et al. [82]. 
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Table 4.8; Dynamic polarizabilities of the OCS molecule". 









EOM-CGSD" 




uj (a.u.) 


TDHF'’ 


CCSD'' 


Cl-like 


Quadratic 


exp. 


.0000 


32.99 


34.69 


(a(0)) 

36.96 


35.71 


34.35*^ 


.0720 


0.90 




1 

CO 


1.09 


0.67" 


.0886 


1.33 




1.70 


1.65 


1.01" 


.0934 


1.56 




1.91 


1.86 


1.21" 


.0995 


1.79 




2.20 


2.17 


1.42" 


.0000 


24.01 


27.55 


A«(0) 

32.50 


29.97 


26.25'^ 


.0720 


1.06 




Ao;(a;) — Ao;(0) 
1.59 


1.48 


1.62" 


.0886 


1.57 




2.42 


2.29 


2.57" 


.0934 


1.85 




2.78 


2.57 


2.83" 


.0995 


2.12 




3.19 


2.97 


3.11" 



“ The experimental bond length for R(OC) = 2.191 a.u. and R(CS) = 2.947 a.u. 
* Present work. 

" G. R. Almset al. [82]. 

^ Extrapolated from G. R. Alms et al. [82]. 
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Table 4.9: Dynamic polarizabilities of the CS 2 molecule". 









EOM-GGSD" 




to (a.u.) 


TDHF'' 


CCSD* 


Cl-like 


Quadratic 


exp. 


.0000 


55.18 


55.03 


(a(0)) 

62.10 


59.15 


55.40“^ 


.0720 


2.50 




{a{uj)) - (q;(0)) 
2.84 


2.89 


3.24" 


.0886 


3.93 




4.67 


4.62 


5.27" 


.0934 


4.42 




6.49 


5.25 


5.94" 


.0995 


5.10 




7.47 


6.05 


6.82" 


.0000 


56.58 


56.50 


Aa(0) 

73.74 


66.51 


58.71^^ 


.0720 


5.03 




Aa(a;) — Ao:(0) 
6.27 


6.24 


6.41" 


.0886 


7.97 




10.34 


10.01 


10.66" 


.0934 


9.00 




12.38 


11.33 


12.28" 


.0995 


10.44 




14.41 


13.19 


14.31" 



“ The equlibrium experimental bond length for R(CS) = 2.944 a.u. 
* Present work. 

" G. R. Almset al. [82]. 

Extrapolated from G. R. Alms et al. [82]. 
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Frequency (a.u.) 

Figure 4.1: Mean dynamic polarizability of CO 2 in the POLl basis set. 



component increases as 0.6, 2.3 and 4.1 a.u in the series CO 2 , OCS, CS 2 , while 
for the anisotropic component this effect is enhanced from 1.0, 2.5 to 10.5 a.u. 
(!). Invariably, inclusion of orbital relaxation effects brings the results closer to 
experiment. Comparing the EOM Cl-like and EOM quadratic models we find 
that the Cl-like version overshoots the experimental results by up to 15 a.u. for 
the anisotropic component of CS 2 . It is clear that for these molecules derivative 
techniques will work better than expectation-value based methods, but most of 
the error corresponds to orbital relaxation, which is not included in either method. 
Interestingly, the TDHF model performs better than the correlated EOM methods 
for these systems. Relaxation effects are, of course, precisely what is accounted 
for in TDHE. This, together with some unknown, but favorable, cancellation of 
errors is likely to account for the relative accuracy of the TDHF results. Of course, 
we also assume that the reference experimental values are reliable, even though in 
many cases these are rather old results. 
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Figure 4.2: Mean dynamic polarizability of CS 2 in the POLl basis set. 

Dispersion effects in EOM-CC, on the other hand, are in reasonable agreement 
with experiment and greatly improve upon the Hartree-Fock dispersion curves. We 
note that the EOM Cl-like model performs a little better that the EOM quadratic 
model compared to experiment, especially for higher frequencies. 

From these comparisons we conclude, that if orbital relaxation effects appear to 
be important, they should be included. The potential importance of orbital relax- 
ation in time-dependent properties is a tricky question. It is certainly possible to 
include orbital relaxation by building the correlated calculation on the underlying 
TDHF solution, just as coupled perturbed Hartree-Fock (CPHF) solutions are es- 
sential in analytical gradients and Hessian evaluations for the total energy. In that 
case orbital relaxation can never be neglected without invalidating the exactness 
of the critical points. The CCSD and higher methods, though, build in the vast 
majority of such relaxation effects via \ 4 >) = e^‘|(/)) [33, 77], and such flexibility is 
usually sufficient for other kinds of properties [80]. However, when orbital relax- 
ation is as large as found here for CS 2 , we have to reconsider whether they shoidd 




+ TDHF 

X EOM-Q 

X EOM-CI 

□ Exp. 
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be included in highly accurate calculations. However, to incorporate TDHF relax- 
ation means also introducing artifactual TDHF (RPA) excitation energies (part 
of the propagator) into the calculation, which is unsatisfactory. Another obvious 
possibility is to include effects from connected triple and higher excitation opera- 
tors in the calculation. Until the full Cl is reached, though, there are still residual 
orbital relaxation effects. Here, apparently that effect is numerically significant. 
Interestingly the dispersion is, in all cases considered, well described by the Cl-like 
model in EOM-CCSD. If we combine this with a finite field CCSD value for the 
static component (or full second order analytical CCSD derivative), we may expect 
to obtain reliable results, even when orbital relaxation effects are important. 

Trans-butadiene 

Table 4.10 shows the results of the dynamic polarizability calculation for the 
trans-butadiene molecule. Polyenes, of which trans-butadiene is a prototypical 
example are of substantial interest in NLO materials [85]. One can notice that 
TDHF values are much larger than the corresponding EOM values, due to the 
overestimation of the xx and yy components of a. The dispersion behavior is 
also very different: the percentage dispersion at a; = 0.0995 a.u. is 10.3% of 
the TDHE level, 7.61% for the EOM Cl-like method and 5.35% for the quadratic 
scheme. The quadratic correction accounts for up to 2-3% of the total correlation 
correction and also slightly changes the dispersion behavior: for the static case the 
difference between the Cl-like and quadratic anisotropy is only 0.75 a.u. whereas 
for the largest computed frequency it amounts to 3.1 a.u. This result suggests that 
the dispersion of the dipole polarizability in conjugated hydrocarbons is sensitive 
to the correlation level. The fact that correlated results have been found to be 
further from experiment than TDHF [5, 73] has led to the argument [86] that for 
butadiene, and indeed linear polyenes in general, TDHF gives a better description 
of the hyperpolarizabilities than correlated methods. If so, the ultimate correlated 
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Table 4.10: Dynamic polarizabilities of the C4H6 molecule". 









EOM-CCSD" 




u) (a.u.) 


TDHF'’ 


CCSD'’ 


Cl-like 


Quadratic 


exp. 


.0000 


52.04 


48.79 


(«(0)) 

49.78 


49.16 


58.31 


.0720 


2.58 




P 

bo 1 


1.16 




.0886 


4.10 




2.92 


1.97 




.0934 


4.63 




3.32 


2.26 




.0995 


5.37 




3.79 


2.63 




.0000 


50.29 


35.24 


Aq;( 0) 

36.26 


35.51 




.0720 


5.01 




Ao;(a;) — Ao;(0) 
2.60 


1.18 




.0886 


8.06 




4.13 


2.20 




.0934 


9.14 




4.66 


2.55 




.0995 


10.67 




5.41 


3.02 





" The experimental equlibrium bond length r(Ci — C 2 ) = 2.508 a.u. , r(C2 — C3) 
= 2.768 a.u., r(C*i — Hi) = 2.028 a.u., r(C2 — i/3) = 2. 035a. u. 

* Present work. 

‘ |87| 



methods must regain the TDHF results. However, it is difficult to see how the 
correlation effects can be neutralized. The extension of the atomic basis set and 
inclusion of vibrational effects must be considered. 

Long-range London dispersion coefficients. 

Table 4.11 lists the Ce coefficients for the collection of atoms and molecules 
considered in this study. We use the same procedure previously discussed [88]. 
The largest errors are for CS2 and OCS, consistent with our observation of the 
quality of the calculated polarizability results of the two respective molecules. 
Inclusion of the quadratic contribution improves the agreement of the calculated 
results with experiment and inclusion of orbital relaxation would presumably lead 
to still further improvements (see previous section). 
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There are several other correlated theoretical studies of Ce coefficients of small 
molecules in the literature [89, 90, 88], the latter using EOM-CCSD. However, 
the only previous theoretical results of Ce coefficients for CS 2 and COS are due 
to Spademan. The first column of Table 11 contains Ce coefficients calculated by 
Spademan [91] at the TDHF level using a 6-31 G basis set with added polarization 
functions. There are no previous theoretical or experimental Ce coefficients for the 
CI 2 molecule available in the literature. 

Conclusions 

The main focus of this chapter is to study the variation of the molecular polar- 
izability with frequency (both real and imaginary) and to calculate the dispersion 
coefficients (Ce) for a series of molecules at a consistent level. For this purpose, 
we use the recently developed FOM-CCSD Cl-like, linear and quadratic approxi- 
mations to calculate frequency dependent polarizabilities of several selected small 
molecules. The FOM-CCSD Cl-like approximation when combined with Sadlej 
polarized basis sets was shown previously to perform well for molecular dynamic 
polarizabilities [21]. In this thesis we also consider orbital relaxation effects on 
calculated static polarizabilities through finite difference CCSD calculations. 

In general the difference between the FOM Cl-like, linear and quadratic mod- 
els is found to be minor, both for the net result and the dispersion. Also the 
results agree reasonably well with experiment. An important exception is the CS 2 
molecule, and to a lesser extent CCS. In these cases we found orbital relaxation 
effects on the calculated static polarizabilities to be very important, implying signif- 
icant potential effects on the dynamic polarizability. Fven in these cases, however, 
the frequency dependence of the polarizabilities is well described at the FOM- 
CCSD level. If orbital relaxation is important a hybrid method like a combination 
of FOM-CCSD for dispersion and finite-field CCSD (or analytical CCSD second 
derivatives) for the static component could offer a pragmatic approach. 
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Table 4.11: Cg dispersion coefficients for homonuclear interaction. 





TDHF* 


EOM-CCSD (CI-like)“ 


EOM -CCSD (Quadratic)" 


exp. 


N2 


71.46 


73.63 


71.84 


73.40" 


CD 


- 


416.3 


405.2 


- 


C 2 H 2 


196.6 


204.3 


196.9 


231.4'^ 


CO 


63.29 


83.51 


80.64 


81.40" 


CO 2 


123.1 


173.6 


166.4 


158.7" 


SCO 


368.5 


453.1 


481.4 


402.2" 


CS 2 


826.0 


1182 


1079 


991.3" 



“ Present work. 

M. A. Spackman [91]. 

Margoliash and Meath [92] and M. A. Spackman [91]. 

Starkschall and Gordon [93]. 

® Fowler et al. [94]. 

To better include orbital relaxation, in particular also the contribution to dis- 
persion, the EOM-CCSD model can be extended with an approximate treatment 
for triples, like EOM-CCSDT-1 or EOM-CCSDT-3 [95]. Such triple contributions 
have been shown to introduce the most important remaining orbital relaxation 
terms [80]. The numerical effect of triple excitations (in cases that relaxation ef- 
fects are not considered very important) on polarizabilities can be assessed from 
other results [77, 96], with the general conclusion that it increases the magnitude 
of the polarizability, but not necessarily the agreement with experiment. 

Dispersion at the TDHF level is often a poor representation of the correlated 
EOM-CC dynamic polarizibility in the same basis set. Dispersion at the EOM- 
CCSD level agrees significantly better with the experimental measurements. 

4.2 First and Second Hyperpolarizability 
4.2.1 Hydrogen Fluoride 

A number of papers have studied the hyperpolarizabilities of the FH molecule 
with correlation since 1979 [97, 98, 99, 100, 101, 102, 103, 85]. A critical analy- 
sis of all the high-level calculations until 1993 [101] concluded that, contrary to 
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several other molecules whose theoretical frequency-dependent hyperpolarizabili- 
ties agreed to within about 10% with gas phase experimental data, the theoretical 
results for FH showed an error of about 20%. Despite the small difference, and 
the very small magnitude of the FH hyperpolarizability that can enhance the per- 
cent error, high-level, correlated ab initio theory should be able to unambiguously 
provide agreement to within 10% for such a simple molecule, or the experimental 
value [104] should be reconsidered. 

In prior theoretical work, the frequency dependence had been estimated from 
time-dependent Hartree-Fock (TDHF) theory [101], or from multiconfigurational 
linear response (MCLR) theory [100]. However, the former is uncorrelated and 
the latter includes mostly nondynamic correlation while much of the magnitude of 
hyperpolarizabilities arises from the large dynamic correlation effects. The latter 
are, perhaps, most accurately and conveniently included via coupled-cluster (CC) 
theory, but except for some Kerr effect studies [47], the prior CC results have not 
included the correlated frequency dependence that is consistent with CC theory. 
In this letter we have generalized the equation-of-motion (EOM-CC) method to 
recursively and analytically compute correlated, frequency dependent hyperpolar- 
izabilities, to any order. This provides an unambiguous correlated measure for the 
dispersion effect. Furthermore, although many different kinds of basis sets have 
been considered which appeared to offer consistent results, there has not been a 
more systematic study based upon a series of n-aug-cc-pVXZ basis [105] to at- 
tempt to farther eliminate this degree of uncertainty in the calculations. Here, we 
present such EOM-CC frequency-dependent results along with several other bases, 
to attempt to find basis set limit results. 

The EOM-CC approach in its CTlike approximation [21, 53, 44] offers a con- 
venient recursive treatment of second [21] and higher-order properties. We sim- 
ply consider the usual perturbation expressions that can be used to define such 
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frequency-dependent properties. ^ Hence, all of the frequency-dependent polariz- 
abilities can be written in the convenient form, 

’ ^2 ) • • •) P (ik^k) (0| (l “f Qii |b) R ( iO\ LO 2 ~ ^n) 

*(h|^2|h)R(-a;2 - * ••• * (h|9j„|h)R(-o;„) (/i|gi„^JO) 

(4.12) 

where P(ifc,o;fc) denotes all the (n-l-1) permutations of the pairs and 

R(<^x) = (h|£'o — Hq H-o; 2 ;|h) ^ (4.13) 

is the resolvent operator for the particular frequency lo^. Also Qi = Qi — (qi), 
where (g,) = (0| (1 -I- A) gj |0) and qi = e~^qie^. Here, A is the deexcitation 
operator introduced previously [42] and qi denotes the i-th component of the dipole 
moment vector. It should be clear from the above that, operationally, the repeated 



^ An alternative definition of second and higher order properties is offered for 
static (a; = 0) cases by derivatives of the energy or the induced dipole moment. As 
is well known, even in the case of an expectation value versus the first derivative, 
the two definitions will not yield exactly the same results unless the Helman- 
Feynmann theorem is satisfied. For first-order properties in CC theory, this equiv- 
alence is ensured by virtue of the CC expectation value < 0|(1 4- A)e“^0e^|O > 
for an operator © where both A and T are determined from stationary equa- 
tions. However, the usual perturbation expansions for higher order properties are 
only strictly equivalent to higher derivatives when the higher-order analogues of 
the Helmann-Feynman theorem are satisfied, and this is not always the case [44]. 
The difference can be estimated from the numerical derivative value (unrelaxed) 
compared to the perturbation theory expression. Here that difference for the t-aug- 
cc-pVTZ basis is about 0.1 a.u. for /?||(0) and 7||(0) so for all practical purposes 
for a single molecule there is no significant difference for this example. For other 
comparisons see [53, 106, 40]. This additional small effect is incorporated in the 
orbital relaxation estimates of this thesis. For the frequency dependent case, the 
perturbation theory definition can be immediately generalized by simply adding 
the appropriate non-zero u! to the resolvent [21] while the usual derivative approach 
is not directly applicable and requires some additional considerations [39, 53]. 
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evaluation of the vectors 

Tg> = {h|«,|h)7’<‘> (4.14) 

and 

R(c^,)7)f (u-J (4.15) 

where 

= R M {h\m (4.16) 

are the critical computational steps for all hyperpolarizabilites in the perturbation 
theory definition. This is actually accomplished by solving the corresponding linear 
equations of which the first is 

R-i(u;,)7)(') = (/r|%|0) (4.17) 

rather than taking the inverse of the resolvent [21]. Hence, there is no truncation 
of the excited states defined over the configurations |h) as is often the case in sum 
over state approximations. Obviously, we can go to any order hyperpolarizability 
with the same program, just as was accomplished in the open ended TDHF method 
presented sometime back [3], making a, /?, 7, 5, e etc., readily accessible for whatever 
optical process, by virtue of using the correct form of R((^j;). 

Using the above, numerical results in several bases are presented in Table 4.12. 
The modest J basis used previously consists of [6s5p5dlf] on F and [5s3pld] on H 
and appears to be the best converged of our prior calculations [101]. The POLl 
basis, due to Sadlej, is defined for polarizabilities [71]. Here we augment it with 
two additional shells of diffuse functions obtained by 0.3 times the most diffuse 
exponent of each in the POLl basis to give [7s5p4d2f] on F and [5s4pld] on H. The 
series of bases cc-pVTZ, aug-cc-pVTZ, d-aug-cc-pVTZ, t-aug-cc-pVTZ, q-aug-cc- 
pVTZ and p-aug-cc-pVTZ i.e. [9s8p7d6f/8s7p6d] with augmentation as considered 
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elsewhere [105] are hoped to provide a converging series. We also consider a t-aug- 
cc-pVQZ basis [8s7p6d5f4g/7s6p5d4f] to offer a different extension. Of course, as is 
usually the case for smaller bases, the FH (5\\{ /?n = 1/5 Ei {Puz + Pzu + Pizi) ) tends 
to be too large, while adding further diffuse augmentation reduces the magnitude, 
which then eventually might grow, gradually, with basis set enhancement. 



First, depending upon the n-augmented cc-pVTZ series, we would expect a 
second-harmonic generation (SHG) result of about -8.8 for a pentuply augmented 
cc-pVTZ basis, a result in line with the others shown. The larger values for \Pshg 
correspond to probably excessive values for \Pzzz\, which is smaller in better basis 
sets. All of these results are at the EOM CCSD level and are still 20% in error com- 
pared to the center of the experimental value [100]. The dispersion correction given 
by EOM CCSD is 9.9%, in excellent agreement with the 9.6% estimate from TDHF 
that was previously used, and the 10% MCLR value. Hence, the prior dispersion 
estimates [101] were entirely justified for this problem , and did not introduce a 
significant error. We see a similar behavior for the direct-current induced second- 
harmonic generation process {^dc-SHc)- Our J basis and extrapolated pentuply 
augmented result is 650 a.u., also 20% in error. 

Besides dispersion and basis, there are several other items previously consid- 
ered [101, 97, 98] to try to resolve the discrepancy between theory and experiment. 
These include the role of triple excitations as measured by static CCSD(T) results, 
vibrational averaging (i.e. to tq) in the ground state [98] , and the pure vibra- 
tional correction [107] that arises from vibronic intermediate states in perturbation 
theory expressions. 

The magnitude of the triple excitations in CCSD(T) static values is -0.69 a.u. 
in the J basis [101], giving a CCSD(T) static value of p\\ equal to -8.72 a.u.. In 
the t-aug-cc-pVTZ basis, CCSD(T) is -8.54 a.u. Using the latter triples effect 
(-0.62 a.u.), one would estimate a p-aug-cc-pVTZ CCSD(T) /5||(0) value of -8.57 



Table 4.12: EOM-CCSD unnormalized static and dynamic {lo = 0.0656 a.u. / 694.3 nm) hyperpolarizabilities of the HF molecule 
(in a.u.) at Rhf=^-7S2S a.u., except where indicated. 
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a.u. (Though the second decimal is certainly not significant in these estimates, we 
carry it temporarily for ease of analysis). Applying the same dispersion obtained 
in the EOM-CCSD calculations to the triples, we would obtain -9.61 a.u. in J 
and -9.48 a.u. for the estimated CCSD(T) SHG process in the p-aug-cc-pVTZ 
basis, assuming that there is no effect of the triple excitations on the dispersion. 
The computed pure vibrational correction is -0.35 a.u. [107] while -0.38 a.u. is 
the currently determined effect in the J or the t-aug-cc-pVTZ basis for evaluating 
EOM-CCSD Pshg using formula (24) from [108] at four different internuclear 
distances^ . Adding these to the EOM-CCSD SHC result gives -10.3 a.u. in J and 
-10.2 in the p-aug-cc-pVTZ basis, slightly within the error bars of the experiment. 
Why, then, is this different from all previous calculations [96] ^ which supported a 
maximum \(3shg\ CCSD(T) result of -9.0, or -9.4, if the pure vibrational correction 
is included? 

The principal difference between the present calculations and the prior ones 
is that in EOM-CC (or the related CCQR treatment of second- and higher-order 
properties [39] ^ ), the orbitals are not allowed to relax to accommodate the per- 
turbation, as they would be in a full derivative calculation. The important role 



^ This effect was previously found to be -0.28 to -0.37 a.u. [101], -0.48 [107], 
with the MCLR value being -0.62 [100]. The four values used in the J basis CCSD 
calculations of the SHC zero point average (using Bishop’s formula [108]) are -8.53, 
-8.85, -9.18 and -9.52 a.u. for HE separation distances of 1.7026, 1.7328, 1.7620 
and 1.7922 a.u. respectively. 

^ Shelton and Rice [102] cited unpublished CCSD(T) results in a large basis give 
/?ll(0)=-7.31 a.u. 7|| (0)=525 a.u. These are nearly identical to the J basis result 
of [96], namely -7.43 a.u. and 525 a.u. 

^ A recent paper [40] reports SHC for HE using CCQR. In the same t-aug- 
ccPVTZ basis and frequency used here, the result is -8.789 a.u. The difference of 
0.089 a.u. is a measure of the effect using the usual perturbation theory expression 
compared to the derivative expression. 
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of orbital relaxation for hyperpolarizabilities has been demonstrated previously 
« 

[80]. This relaxation effect can be estimated by a finite field static calculation. 
This effect can be quite large as shown in Table 4.13. In the J Cartesian basis, 
the reduction in /3||(0) is from -8.03 to -6.89(1). With the triples (whose relaxed 
value is -0.54 a.u.), we see how the best previous CCSD(T) /?n(0) value of -7.4 a.u. 
arises. If we assume the same relaxation for the SHG process, our J basis result 
would be -7.71 which is close to the estimated CCSD J basis result before [101], 
and is more than 20% in error. Adding the above triples with estimated disper- 
sion and vibrational corrections (i.e. -1.27 a.u.), our J basis SHG result would be 
-9.0 a.u., consistent with the best previous estimate. Doing the same exercise for 
the d-aug-cc-pVTZ, we obtain a relaxation effect of -1-0.52 in the spherical basis, 
giving a result of -9.3 a.u. If we also add the change due to the extra Cartesian 
components, assuming that in the absence of linear dependencies (none are found) 
the more functions, the better, we get a final extrapolated result of -9.1 a.u., quite 
close to that in the J basis. In the t-aug-cc-pVTZ case the Cartesian and relax- 
ation effect is 4-0.55 a.u. (two linearly dependent functions are removed in the 
calculation). That makes the result become -9.4 a.u., still outside the error bars, 
but identical to the previous \(^shg\ maximum of -9.4 a.u. The extrapolation to 
p-aug-cc-pVTZ gives -8.8 4- (-1.3) = -10.1 a.u., while geometric extrapolation of 
the relaxation/Cartesian effect in the p-aug-cc-pVTZ basis suggests 4-0.3 a.u. to 
return to -9.8 a.u. still slightly outside the error bars, but better than the previous 
-9.4 a.u. estimate. Since we certainly believe that barring pathological situations, 
’’relaxed” orbitals, are better than unrelaxed ones, and particularly for a property 



like hyperpolarizabilities that the more functions the better, which favors Cartesian 
bases over spherical ones; the best J, d-aug-cc-pVTZ and t-aug-cc-pVTZ results 
are consistent with all prior calculations, while the more extensive extrapolated 
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Table 4.13: CCSD and ( CCSD(T) ) orbital relaxation and basis set effects (in 
a.u.) in the j3 and 7 static hyperpolarizabilities of HF 



Relaxed orbitals Unrelaxed orbitals 



Basis Spherical Cartesian Spherical Cartesian 



• 

J 

p\\ 


-7.50 


-6.89 (-7.60) 


-8.06 


-8.03 (-8.72) 


7|| 


- 


478 (526) 


- 


574 (631) 


d-aug-cc-pVTZ 


P\\ 


-7.24 


-7.07 


-7.76 


-7.63 


7|| 


(538) 


- 


(591) 


- 


t-aug-cc-pVTZ 




-7.51 


-7.37 


-7.92 


-7.84 


7|| 


530 (582) 


539 


558 (613) 


571 


q- aug- cc- p VT Z 


/^ll 


-7.60 (-7.87) 


- 


-7.95 (8.19) 


- 



p-aug-cc-pVTZ result suggests some improvement to -9.8 a.u. Even the use of the 
very large t-aug-cc-pVQZ basis does not change this result. 

We can make the same type of estimate for the experimentally known dc-SHG 
process for 7. The previous theoretical results (with exclusion of the pure vibra- 
tional correction) were estimated to show a value of 660±60 a.u. [101] compared 
to 840±120, experimentally.' 

The best of the new EOM-CCSD results, including the extrapolated p-aug- 
cc-pVTZ and previously used J basis are 650 a.u. The EOM-CCSD dispersion 
effect for dc-SHG is 12.7% in the t-aug-cc-pVTZ basis compared to TDHF value 
of 12.5% and much higher MCLR value of 19% [100]. For 'fdc-SHGy we can es- 
timate the contributions from triples corrections, vibrational averaging and the 
pure vibrational part from our current work and from [107]. In the J basis triples 
corrections amount to 57 a.u. plus 8 a.u. for dispersion (65 a.u.). Adding the pure 
vibrational correction of -5.5 a.u. from [107], and 3 a.u. computed by changing 
to ro, without orbital relaxation, we get 62 a.u. added to the 'jdc-SHG EOM-GCSD 



result, to give 712 a.u. as the estimated CCSD(T) 'ydc-SHG in the J basis, just 
outside the experimental error bars. If we include the very large static CCSD(T) 
relaxation (-105 a.u.) in the J basis, we get 607 a.u., well outside the error bars. In 
the d-aug-cc-pVTZ basis the relaxation is reduced to 53 a.u. In the t-aug-cc-pVTZ 
basis, the effect of triples is 55 a.u. plus dispersion (62 a.u.), but the relaxation 
effect is only -31 a.u., giving 663 a.u. The Cartesian components would appear 
to add 9 a.u. plus - 5.5 a.u. for the pure vibrations to give 667 a.u. Finally, 
the estimated p-aug-cc-pVTZ CCSD(T) dc-SHG (extrapolation of relaxation ef- 
fect suggests 10 a.u.) gives a nearly identical 666 a.u., again well outside the error 
bars. Despite the best efforts from 18 years of quantum chemistry, it would appear 
that there is still some disagreement between theory and experiment for (5s hg ^iid 
'Ydc-SHG for the FH molecule. The best current theoretical results for 'Ydc-SHG and 
fisHG do not quite fall within the rather generous error bars. If orbital relaxation 
could be neglected, agreement would be achieved even in modest basis sets. In a 
complete basis, of course, there is no relaxation effect, but other than that, there 
are two arguments for dismissing orbital relaxation, one practical and one formal. 
Practically, when using an RHF reference in a correlated calculation, at large bond 
lengths, relaxed RHF orbitals have the wrong behavior, since RHF goes to an ionic 
separated atom limit, which the correlated method is attempting to undo. How- 
ever at equilibrium, even for a fairly ionic molecule like HF, one would not believe 
that unrelaxed orbitals are closer to the truth than relaxed orbitals. The formal 
objection is that inclusion of Hartree-Fock orbital relaxation introduces artificial 
excitation energy poles in the dynamic polarizability, in addition to the correct 
ones [109] 

However, it is absolutely necessary to use relaxed orbitals in gradient and hes- 
sian calculations, as the critical point would not be well defined otherwise. Why 
should it be different for other kinds of properties? The orbital relaxation terms 
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introduced by adjusting the orbitals in the presence of the perturbation would 
correspond to an infinite sum of diagrams if they were introduced via evaluation 
in terms of the unrelaxed orbitals [110, 80]. CCSD easily sums many of these 
diagrams [80], (CCSDT sums more) but not all. For many properties, the infinite 
sums inherent in CC makes orbital choice largely irrelevant [111], but apparently 
not for the hyperpolarizabilities of FH, at least at the CCSD(T) level, as this subtle 
effect accounts for the current and previously observed differences between theory 
and experiment. In a complete basis, relaxation is irrelevant, but we do not seem 
to approach that in the current calculations, even in the p-aug-cc-pVTZ basis. For 
the more common applications to polyatomic molecules, it is even less likely, and 
it would appear orbital relaxation should not be neglected in rigorously computed 
answers. Here, our best final, ’’relaxed” results are -9.8 a.u. for Pshg and 670 
a.u. for Jcic-SHG- Considering the estimates made, a ± 5% error would probably 
be justified in each (-9.8 ± 0.5 a.u. and 670 ± 40 a.u.), causing the experimental 
and theoretical error bars to overlap. However, the center of the experiment and 
theoretical calculation remains about 11% apart for Pshg and 19% for 'fdc-SHG- 
If orbital relaxation could be neglected, our best results would be -10.2 ± 0.5 a.u. 
and 700 ± 40 a.u. compared to -11 ± 1 and 840 ± 120, experimentally. 

4.2.2 fran^-Butadiene 

Trans-butadiene has served for many years as the prototype system for polyenes 
which have attractive intrinsic nonlinear optical properties (NLO). Clearly, if quan- 
tum chemistry is to play a role in the design of NLO materials, it must be able to 
correctly account for butadiene’s hyperpolarizability. 

However, it has been pointed out in several papers [5, 47, 85] that trans- 
butadiene’s theoretically determined hyperpolarizability behaves differently com- 
pared to several other molecules [101, 85]. Whereas for N 2 , CO, CO 2 , H 2 , HF, 
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HCl, NH3, H2O and others, correlation effects evaluated with MBPT or CC the- 
ory provides substantial improvement compared to experiment [101, 85], as they 
should, for butadiene correlation appears to hurt the agreement, as the uncorre- 
lated TDHF value for 7|| for butadiene (and to a lesser extent ethylene) appears 
to give very good agreement with experiment while MBPT(2) and CCSD appears 
to change this value dramatically [5, 47]. Furthermore, recently, Norman et al. 
[112] made the statement with reference to TDHF that ”... none of the more so- 
phisticated correlating methods are able to match this accuracy...,” implying that 
TDHF applied to dynamic hyperpolarizabilities of frans-butadiene and other lin- 
ear polyenes provide accuracy exceeding that of highly-correlated coupled-cluster 
methods. 

The history of quantum chemistry is checkered by many examples of excellent, 
but accidental agreement between theory and experiment, largely based on fortu- 
itous cancellation of errors. If TDHF ’’gets the right answer for the right reason,” 
then all correlation effects must effectively vanish, which seems highly unlikely. To 
answer this question, a more detailed study is necessary, taking into consideration 
basis set effects, vibrational corrections, higher excitation contributions and molec- 
ular geometry. Although it has been possible to do MBPT and coupled-cluster 
calculations of static hyperpolarizabilities for a long time [97, 5, 101], a consistent 
approach to dynamic processes has not been possible until recently [47, 48, 40]. 
We can now calculate frequency dependent hyperpolarizabilities for any order and 
any process, recursively, using the EOM-CCSD method [37, 68, 20, 21]. Hence, we 
are now in a position to offer more definitive correlated results. 

The average hyperpolarizability of trans-butadiene in C2h symmetry for the 
dc-SHG process is defined as 



7|| 
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The actual working equation for a single component of the dc-SHG process is quite 
simple when based upon the standard perturbation expression [21, 48]: 

72222 (-2o;;a;,a;,0) = P {oj) ](0| (1 + A) (-2a;) (-u;) (0) 9 ^| 0 ) - 

(0| (1 + A)q^R{-2uj)qz\0){0\ (1 + A) R{-u) qzR{uj) q^\0) - 
(0| (1 + A) q,R (-2a;) gA|0)(0| (1 + A) q,R (-a;) R (co) q,\0)] 

(4.19) 

where the permutation operator P(o;) permutes the frequencies (—2a;, a;, a;, 0) with 
their associated spatial subscripts (here: z,z,z,z). Similarly, we can readily eval- 
uate the THG, IDRI and EOKE processes. In the above expression Ro{nu>) is a 
frequency-dependent resolvent: 

R (±no;) = |h)(h|E'cc — Hq± no;|h)“^(h| (4.20) 

In our current calculations we employ three sets of basis functions: that previ- 
ously used in irans-butadiene calculations [113, 5, 47] 6-3lG-|-(P,D) basis [3s3pld/2s] 
obtained by augmenting the standard 6-31G set by one p and one d gaussian with 
the exponent C=0.05, the standard POLl set [71] [5s3p2d/3s2p] and finally, the 
doubly augmented POLl set [5s4p3dlf/5s4pld], denoted POLl-f-f [48] here. The 
geometry plays an important role in these calculations, especially the TDHF ones. 
We choose the three sets of coordinates in Table 4.14: that obtained by the SCF 
optimization in 6-31G basis (geometry I) previously used, that obtained by a GCSD 
optimization in the 6-3lG-|-(P,D) basis (geometry II), and the ’’experimental” ge- 
ometry from Haugen and Traetteberg [114]. 
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Table 4.14: Geometry of the ^ran^-butadiene molecule at various levels of opti- 
mization 



Variable 


SCF (I) 


CCSD (II) 


exp“ (HI) 


C 2 -C 3 


1.465A 


1.481A 


1.467A 


C 1 -C 2 = C 3 -C 4 


1.327A 


1.361A 


I. 343 A 


C 1 -C 2 -C 3 


124.3° 


124.0° 


122 . 8 ° 



“ Ref. [114] 

In the small basis set originally used, 6-3lG-l-(P,D) [5, 47] the static EOM- 
CC averaged hyperpolarizability is 19.61 (here we assume the units of 10^ a.u.), 
whereas the corresponding SCF result is 13.62. Our previous studies [48] show a 
large effect of including the Cartesian component of the d shell in smaller basis 
sets. Introducing this component the EOM-CC and SCF hyperpolarizabilities are 
20.17 and 14.81, respectively. 

The essence of the controversy is that for dc-SHG the TDHF results in the 
6-31G+(P,D) Cartesian basis, at o;=0.0430 a.u. and w=0.0656 a.u. are 18.2 and 
24.9, which are in reasonable agreement with experiment (20.2 and 27.7 respec- 
tively). The above static MBPT(2) correlated results in this basis augmented by 
the percentage TDHF dispersion to estimate the correlated dc-SHG results, were 
substantially in error at 25.3 and 34.8 respectively [5]. Of course, there are many 
potential sources of error in such a calculation, particularly the dispersion esti- 
mates, and a different estimate of the dispersion [47], based upon the relation of 
Bishop [115], plus CCSD instead of MBPT(2), reduced the correlated value to 24.3 
and 29.7, in better agreement with experiment. However, neither of these results 
is satisfactory, as the proper result comes from a correlated, frequency dependent 
calculation. In order to obtain a more definitive result, we consider much larger 
basis sets and geometry effects, as well as the correlated frequency dependence. 

Table 4.15 presents 6-3lG-l-(P,D) TDHF and EOM-CCSD frequency dependent 
hyperpolarizabilities in the three different basis sets at geometry I. Since the POLl 
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Table 4.15: Static and dynamic hyperpolarizabilities of trans-butadiene (in 10^ 
a.u.) in geometry I. 



In parentheses, values obtained with all 6 Cartesian components of the d shell. 







THG 


dc-SHG 


IDRI 


EOKE 


0.0000 


TDHF 

6-3lG+(P,D) 


13.62(14.81) 


13.62(14.81) 


13.62(14.81) 


13.62(14.81) 




POLl 


13.92 


13.92 


13.92 


13.92 




POL1++ 


13.70 


13.70 


13.70 


13.70 


0.0430 


6-3lG+(P,D) 


20.78(22.96) 


16.53(18.2) 


15.45(16.89) 


14.49(15.80) 




POLl 


22.33 


17.29 


16.02 


14.91 




POL1++ 


22.37 


17.14 


15.84 


14.71 


0.0656 


exp 

6-3lG+(P,D) 


55.21(62.58) 


20.18“ 

22.46(24.93) 


18.55(20.47) 


15.80(17.28) 




POLl 


6149 


24-44 


19.73 


16. 42 




POL1++ 


74.50 


24.60 


19.67 


16.25 


0.0882 


exp 

6-3lG+(P,D) 


35.87 


27.70“ 

40.92 


25.03 


18.02 




POLl 


140.14 


49.24 


28.14 


19.03 




POL1++ 


40.72 


51.67 


28.59 


18.94 


0.0000 


EOM-CC 

6-3lG+(P,D) 


19.61 


19.61 


19.61 


19.61 




POLl 


17.32 


17.32 


17.32 


17.32 




POL1++ 


18.17 


18.17 


18.17 


18.17 


0.0430 


6-31G+(P,D) 


27.00 


22.59 


22.00 


20.74 




POLl 


39.27 


20.57 


19.55 


19.11 




POL1++ 


36.74 


22.45 


20.56 


19.29 


0.0656 


exp 

6-3lG+(P,D) 


40.81 


20.18“ 

29.19 


26.22 


22.51 




POLl 


36.46 


26.63 


23.32 


20. 77 




POL1++ 


39.31 


27.76 


24.57 


20.99 


0.0882 


exp 

6-3lG+(P,D) 


93.63 


27.70“ 

43.68 


33.24 


26.27 




POLl 


101.31 


36.82 


29.48 


23.04 




POLl-h+ 


104.40 


45.50 


31.35 


24.76 



“ Ref. Shelton and Rice [116], Ward and Elliott [117] 



82 



Table 4.16: Static and dynamic hyperpolarizabilities of trau5-butadiene (in 10^ 
a.u.) in geometry II. 

Upper row: TDHF values. Lower row: EOM-CCSD values 



process 


6-3lG+(P,D) 


POLl 


POL1++ 


exp°' 


static, a; = 0 


16.42 

18 . 43 


15.46 

17.19 


15.29 

18.06 


- 


dc-SHG, u = 0.0430 


20.41 

22.04 


19.54 

20.95 


19.47 

21.18 


20.18±0.11 


dc-SHG, uo — 0.0656 


29.01 

29.73 


28.66 

28.61 


29.04 

29.56 


27.70±1.60 



“ Ref. Shelton and Rice [116], Ward and Elliott [117] 

and POLH-+ sets do not show such large differences, we limit ourselves to the 
spherical components only. Eor TDHE the 6-3lG+(P,D) results always lie below 
the experimental values and it appears that the POL1++ basis is converged since 
the difference between POLl and POL1++ is always quite small. On the other 
hand, the EOM values tend to oscillate with a minimum at the POLl values. 

In order to assess the geometry effect the dc-SHG values given in Table 4.15 
for geometry I may be compared with those for geometries II and III given in 
Tables 4.16 and 4.17 respectively. (The TDHE values of 7(0) in the 6-3lG+(P,D) 
basis have a quite large range: from 14.81 at the SCF optimized geometry to 
16.42 at the CCSD geometry.) The EOM results vary very little with geometry 
being 18.9 ± 0.5. As the basis size increases from POLl to POL1++, the TDHF 
static hyperpolarizabilities appear to be well converged, while EOM results increase 
by almost 5% in all geometries. Our EOM-CCSD values of 7(0) in the POLl 
basis are close to those obtained by finite differentiation of the CCSD energy 
by Kirtman et al. [118]. Our results of 17.32, 17.19 and 17.24 in geometries 
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Table 4.17: Second hyperpolarizabilities of trans-butadiene at the experimental 
geometry (III). 

Upper row: TDHF values. Lower row: EOM-CCSD values 



process 


6-31G±(P,D) 


POLl 


POLl±± 


exp°' 


static, a; = 0 


16.09 

18.60 


14.62 

17.24 


14.52 

18.13 


- 


dc-SHG, LO = 0.0430 


20.00 

22.09 


18.33 

20.26 


18.28 

21.27 


20.18±0.11 


dc-SHG, LO = 0.0656 


28.42 

30.23 


26.42 

27.96 


27.54 

29.63 


27.70T1.60 



“ Ref. Shelton and Rice [116], Ward and Elliott [117] 

I, II and III, respectively, are close to the value of 18.0 obtained therein. Hence, 
these numbers show little effect of the non-Hellmann-Feynman terms, that are 
due to orbital and additional CC amplitude response [44], that are excluded here, 
but would be included in the numerical derivatives as discussed elsewhere [48]. 
Table 4.15 contains the values for four fourth-order processes: third-harmonic 
generation (THG) 7(— Sea; ta, ca, ca), dc induced second-harmonic generation (dc- 
SHG) 7(— 2a;; ca, lo, 0), intensity dependent refractive index (IDRI) 7(— a;; a;, —u>, co) 
and the electro-optic kerr effect (EOKE) 7(— o;;a;, 0, 0). The dc induced optical 
rectification (dcOR) 7(0; lu, — o;, 0) after averaging is the same as EOKE so it will 
not be considered separately. 

For 1,3- trans-butadiene the experimental values have been determined only for 
the dc-SHG process [116], [117] by the dc electric-field induced optical SHG exper- 
iment. The effect at the ruby-laser frequency (0.0656 a.u.) is 27.7 ±1.6 while at 
0.0430 a.u. it is 20.2 ± 0.1. In the SGF optimized geometry (Table 4.15) the TDHF 
values are consistently ca. 3 units too small compared to experiment, while EOM 
values are larger by about 1 to 3 units. At the CGSD and experimental geometry 
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the EOM values for dc-SHG are systematically larger than experiment whereas 
the TDHF values may be either larger or smaller. Moving from the 6-31G+(P,D) 
basis to POLl generally improves the EOM results. Adding extra functions to the 
POLl basis increases the EOM-GCSD results by 1 to 2 units. For trans-butadiene, 
the only previous correlated frequency dependent hyperpolarizability calculation 
was evaluated by finite field differentiation of the frequency dependent linear po- 
larizability [47] to give values for the Electro-Optic Kerr Effect (EOKE): they are 
20.7, 21.9 and 23.7 at frequencies o;=0.0, w=0.0430 and o;=0.0656 respectively in 
the 6-31-t-(P,D) basis. Our current results: 19.6, 20.7 and 22.5 are within 6% of the 
previous values and again show little importance of the non-Hellmann-Feynman 
terms that are neglected in the analytical EOM-CC Cl-like expression [48]. 

Figures 4. 3,4. 4 and 4.5 show the frequency dependence of the relevant fourth- 
order processes of trans-butadiene in all three basis sets in the 6-3lG-t-(P,D) op- 
timized geometry. The TDHF dispersion is too high (27% in the POL+-I- ba- 
sis compared to 17% for EOM at t<;=0.0430 a.u. and 90% compared to 64% at 
o;=0.0656 a.u.), which would also cause an estimated percentage dispersion effect 
based upon TDHF superimposed on static correlated results to overshoot. That 
is primarily why the original MBPT(2) estimate of dc-SHG for butadiene was too 
large at 38.4 [5], or the CCSD result at 34.8 [47]. The alternative estimate of 
dispersion provided by Bishop’s formula [115] based upon EOM-GCSD results for 
OKE [47] led to a 29.7 value. The formula is 

7|| ^1) ^ 3 ) = 7 (0) (1 T T + ...j , 

= <^^ + -f -h cjg. (4-21) 

For small frequencies we can set B=0 and evaluate A in a straightforward manner. 
For the EOM method, at ^<-’=0.0656 a.u., the A coefficient ranges from 18 to 22 
depending on the process and basis set. The TDHF values are on the other hand 
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scattered widely, especially for multiphoton processes like THG. The TDHF A 
value for EOKE is 22 but for THG in the POL1++ basis it jumps to 86, making 
any predictions based on Eq. (4.21) with TDHF completely unreliable. A similar 
trend can be observed at any frequency studied: at ci;=0.0430 a.u. TDHF values of 
A range from 17 to 29 while EOM is more consistent placing them between 14 and 
17. At the 0.0430 a.u. field strength A varies between 14 and 17, that is 18%, while 
TDHF varies between 17 and 29, giving a 70% error. This means that if we try to 
use the A coefficient extracted from the TDHF EOKE value (smallest A) and use 
Bishop’s formula to calculate the third harmonic generation (largest value of A), 
we get an answer with such an error. The 18% error of EOM-GC is also large, but 
we would expect to reduce it by including the quartic term Bcof. For the higher 
frequency, cj=0.0656 a.u. the inferior behavior of TDHF is even more evident: 
while the error in estimating the computed frequencies in the EOM-CCSD is 22%, 
TDHF variations of 22 to 86 units give 291% error for TDHF. Clearly, TDHF fails 
to deliver the accuracy expected from a modern ab initio method. 

To assess two additional important effects; relaxation and triple excitations, we 
performed a series of finite field differentiations of the dipole moment (Table 4.18). 
In the 6-31G+(P,D) basis at geometry I orbital relaxation in CCSD decreases the 
7(0) value by about 7.9% to 18.2. Calculating the CCSD(T) hyperpolarizability 
with relaxed orbitals gives 19.2, which is 2.1% smaller than the EOM-CCSD un- 
relaxed figure. For dynamic processes values of these effects only offer estimates. 
The triples correction was scaled by the percentage dispersion of the EOM-CCSD 
result to account for frequency dependence, while the orbital relaxation was left 
constant for all frequencies and geometries. 

When the effect of relaxation and triples is included, along with the pure vi- 
brational corrections [118], one gets the POL1-I--I- values for dc-SHG (we exclude 
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Table 4.18: Adjusted values of dynamic hyperpolarizabilities (dc-SHG) for trans- 
butadiene in the POLH-+ basis. 





EOM-CCSD 


vibrations^' 


triples" 


relaxation^ 


total‘s 


u = 0.0430 
u) = 0.0656 


22.45 

27.76 


-0.50 

-0.40 


geom I 
0.84 
1.09 


-1.55 

-1.55 


21.24 

26.90 


(j = 0.0430 
uj = 0.0656 


21.18 

29.56 


-0.50 

-0.40 


geom II 
0.84 
1.09 


-1.55 

-1.55 


19.97^ 

28.70^ 


u = 0.0430 
UJ = 0.0656 


21.27 

29.63 


-0.50 

-0.40 


geom III 
1.84 
1.09 


-1.55 

-1.55 


20.06^ 

28 . 77 " 



“ Kirtman et al. [118]. 

Estimated from the finite derivatives of the dipole moment in the 6-3lG-|-(P,D) 
basis at geometry I, by the difference between n for CCSD(T) and CGSD aug- 
mented by the EOM-CCSD percentage dispersion correction, 15.1% at 0.043 a.u. 
and 48.8% at 0.0656 a.u. 

Estimated from the finite derivatives of the dipole moment in the 6-3lG-|-(P,D) 
basis at geometry I, at CGSD. 

These should be compared to experimental values: 20.18± 0.11 at a;=0.0430 
a.u. and 27.70± 1.60 at cu=0.0656 a.u. 

^ It has been pointed out to us by the authors of Ref. 17 that their estimates 
for the correlated values (using a slightly different geometry), are 21.2 and 29.8 
at o;=0.0430 and at u =0.0656, respectively. These results are in good agreement 
with ours. 
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Table 4.19: Excitation energies of ^rans-butadiene at the experimental geometry 

Upper row: TDHF values 
Lower row: EOM-CC values 



mode 6-31G(P,D) POLl POL1++ CCSD“ CCSET-S" ^ 



2Ag 


7.44 

6.93 


7.23 

7.03 


7.07 

6.98 


7.23 


6.85 


7.3 


1A„ 


6.85 

6.41 


6.54 

6.59 


6.47 

6.57 


6.61 


6.59 


6.66 


IB^ 


6.65 

6.22 


6.21 

6.30 


6.16 

6.27 


6.40 


6.37 


6.27 


1B„ 


6.07 

6.62 


5.94 

6.37 


5.92 

6.34 


6.42 


6.35 


5.92 



“ In an augmented ANO basis [119] 
* In an augmented ANO basis [120] 
^ Diarmid [121] 



here any zero point vibrational correction) shown in Table 4.18. The excellent 
agreement of these values with experiment is evident. TDHF hyperpolarizabilities 
already contain the relaxation effect, while subtracting the pure vibrational con- 
tribution will give 16.6, 19.0, 17.8 and 24.2, 28.6 and 27.1 at the three geometries. 
As an indirect measure of accuracy for hyperpolarizabilities, we can consider some 
of the low-lying excitation energies provided by different methods. Table 4.19 con- 
tains the results: the upper (roman font) row contains uncorrelated TDHF values, 
the lower, italicized, corresponding to the EOM-CC energies. At the experimental 
geometry the TDHF seems to offer better agreement with the experiment than 
EOM-CCSD, especially in the smallest basis set, the 6-3lC-l-(P,D). Even the no- 
toriously difficult Rydberg- valence mixed state [122] has an error of only 0.15 
eV, compared to 0.70 eV error on the EOM-CCSD level. Only one CC result ob- 
tains less than 6.0 eV [the two-determinant CCSD (TD-CCSD) result is 5.88 eV 
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[123] at a different geometry], but it is unlikely that TDHF gets this result for the 
right reason. Watts et al. [119] have calculated partial triple excitation contribu- 
tions to the excitation energies of frans-butadiene in an ANO basis set. These are 
shown in Table 4.19 for the iterative EOM-CCSDT-3 method (J.D. Watts, private 
communication) to help to assess any significant role for triple excitations. The 
best results ANO/CCSDT-3 and POLl-f-f/CCSD are quite consistent. 

In this chapter we have presented EOM-CCSD correlated, frequency dependent 
results for four non-linear optical processes and compared those to the correspond- 
ing TDHF results. We show dispersion curves up to 0.0882 a.u. Furthermore, we 
estimate the numerical effects of triple excitations while assuming constant orbital 
relaxation (this effect would be expected to differ at different frequencies and would 
vanish at the basis set limit), so we have to avoid claiming too high an accuracy 
in these calculations. However, the estimates affect the CCSD results by no more 
than ± 2 unit, which seems a reasonable error for additional basis set, geometry 
and correlation effects, too. Hence, at any estimated equilibrium geometry, the 
adjusted EOM-CCSD dc-SHG results average 20.4±2 at o;=0.0430 a.u. compared 
to experiment, 20.18± 0.11 and 28.1T2 at u;=0.0656 a.u. compared to 27.70T 
1.60. The corresponding TDHF results average 17.8 and 26.6 respectively, but are 
more accurate at the CCSD and experimental geometries. 

Although the final averaged values are not too different, the TDHF being some- 
what low and the EOM-CC being slightly high, compared to the rather old exper- 
imental value, the consistency of the EOM-CC with basis and geometry, and with 
it offering more consistent A factors in Eq. 4.21, attest to the greater reliability of 
the correlated result. As a ’’Pauling point,” though, the TDHF results might have 
some merit for large, linear polyenes. 
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in the 6-31G+(p,d) basis 




Frequency (a.u.) 



Figure 4.3: Average dynamic hyperpolarizability of ^mn 5 -butadiene in the 6- 
31G+(P,D) basis set. 
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Figure 4.4: Average dynamic hyperpolarizability of trans-butadiene in the POLl 
basis set. 
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Y in the POL1++ basis 

‘av 




Frequency (a.u.) 



Figure 4.5: Average dynamic hyperpolarizability of trans-butadiene in the 

POLl-h-f basis set. 



CHAPTER 5 

CONCLUSIONS AND FUTURE WORK 
In this dissertation we derived and applied two new methods for advanced, 
highly correlated treatment of higher-order electrooptic properties of molecules. 
Both methods are purely analytical and can be used for both static and dynamic 
(time dependent) properties. With a set of justified approximations, the resulting 
formulas are recursive to any order in external perturbation. The first method is 
based on the Rayleigh-Schrodinger perturbation theory and utilizes the coupled- 
cluster correlated ground state and resolvent operator working in the space of 
singly and doubly occupied configurations. We refer to it as RS-EOM-CCSD. 
The advantages of this method, besides the obvious advantages of the standard 
EOM-CC methods are: ease of implementation, recursiveness and readiness for 
application to other than electrooptic processes. The equations for this method 
derived here are in principle exact and we introduce all the approximations based 
on their numerical significance and difficulty of implementation. 

The current development of quantum chemical methods is diverging in two op- 
posite directions. Applicability and accuracy are very often two mutually exclusive 
words in quantum chemistry: the more accurate a method is, the larger amount 
of computational power is needed to tackle a problem limiting potential usefulness 
of such a method for meaningful chemically interesting applications. 

Similarly, the methods presented in this work can be extended in both direc- 
tions. To achieve results which are closer to the experiment, the next step involves 
including higher excitations in the coupled cluster wavefunction. Adding approxi- 
mate triple excitations in the EOM-CC method [95] either as a non-iterative vari- 
ant [EOM-CCSD(T)] or iterative (EOM-CCSDT-3) proved very advantageous for 
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improving the overall accuracy of the excitation energies calculations, especially 
in the case of double excitations. Extension of these algorithms to calculations 
of the polarizabilities and hyperpolarizabilities is straightforward and is expected 
to improve the numerical accuracy significantly. Besides decreasing the numerical 
error, inclusion of approximate triple excitations has been shown [124] to overcome 
an inherent single-reference burden of the EOM-CC method. 

In order to make the presented methods more suitable for tackling real life 
chemical problems, the RS-EOM-CC and BW-EOM-CC can be easily extended 
to calculate magneto-optic properties, like Cotton-Mutton effect (hypermagnetiz- 
ability) [125, 126] or Earaday effect also called circular birefringence. The working 
equations presented in this thesis need to be extended to include magnetic dipole 
moment integrals and GIAO orbitals. Linear magnetic properties (NMR shielding 
and spin-spin coupling constants) have been already implemented by Perera et al 
[127, 128] in the most recent version of the ACESII program package [74] thus 
making the extension trivial. 

As a large portion of the nonlinear optics involves two- or three- dimensional 
extended systems (polymers, surfaces and crystals) an obvious extension of the 
proposed method is to include the translational symmetry and boundary condi- 
tions. The first correlated and fully analytical way of treating one-dimensional 
infinite systems has been recently implemented in our group by Sun and Bartlett 
[129] for the MBPT second-order wavefunction. Both total energies and analytical 
gradients have been reported and the work on coupled-cluster formalism is on the 
way. 

In order to speed up the calculations of the reduced resolvent, which consti- 
tute the major bottleneck in evaluating the recursive EOM-CC equations for (hy- 
per)polarizabilities one can implement the partitioning scheme for the H Hamil- 
tonian. Such a technique has been formulated and implemented by Cwaltney et 
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al [69] and Nooijen et al [70] for excitation energies and indirect nuclear spin-spin 
coupling constants respectively. 

BW-EOM-CC will also be used to calculate the second- and third-order non- 
linear optical responses in excited states of molecules. This new and very promis- 
ing subject has just recently achieved a great deal of attention after experiments 
showed that the nonresonant hyperpolarizabilities originating from excited states 
can be by several orders of magnitude larger than those originating from the 
ground state. Several groups reported both experimental [130] and theoretical 
[131, 132] results for molecules like diphenylhexatriene, fran^-octatetraene and 
smaller polyenes. When pumped with a laser beam at a frequency of electronic 
resonance, the molecule can maintain a large excited state population for times 
long enough to perform nonlinear optical measurements. Because this process is 
nonresonant and involves virtual electronic transitions, one obtains large ultra- 
fast nonlinear optical responses with minimal background absorption. Detailed 
theoretical research in this field could help develop new materials with nonlinear 
susceptibilities many times larger than any material used nowadays. 



APPENDIX A: 

CANCELLATION OF THE DISCONNECTED TERMS IN THE BW-CCSD 

METHOD 

In the CCSD theory with single and double excitations the T amplitude equa- 
tions have the form: 



{D“‘|W - £|V>ccd) = 0 ^ 

(A?l^|0> + + (Atl" (^r|) |0) = E{A“*IV'ccd> 

{D«‘|//| 0 ) + Y. + {D'!^\H\Xqe) = 

k>lc>d 

Y (kl\\cd)titf, = E - E „ Ftfi , ( 1 ) 

k>lc>d 






($“|C (l + Ti + T 2 + l/2!Tf + TiT 2 + l/3!Tf) |$o) , 




{Eo - Ei) ($g’|T 2 |$o) = = 



(<I>g’|C (1 + Ti + T 2 + l/2!Tf + TiT 2 + l/3!Tf + T 1 T 2 + 1/2!T2 + 1/4\T^) |$o) • 




For the singly excited T amplitudes the above equation can be represented as 
a sum of the connected and disconnected elements, as we do not work within the 
linked cluster approximation anymore. 

Thus we have 



(Eo-Ei) mTil^o) 
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where the disconnected part can be expressed in the form: 



- mV (i/ 2T2 + TiTs + l/3!Tf) |$o)dc 



( 5 ) 



From the above we can infer that 



($“|Fe^‘+^^|$o)DC = ^-AF;o, 



( 6 ) 



with being the amplitude coefficients. After substituting the above into Eq. (4) 
we get 

{Eo - Ei) + t^AEo . (7) 

Now, when the t amplitudes are completely converged, on the left hand side 
of Eq. (7) and in the righmost term of this equation are the same and we can write 



{EoAE - Ei) = 



($«|Ve^i+T,|^o)c 




The above equation contains only connected diagrams and is therefore linked and 
size consistent. The denominator {EqAE - Ei) is exactly the BWPT denominator 
(3.101) containing the exact energy. 

Similarily we can prove the size consistency of the doubles equations by express- 
ing the amplitude equations (1) as a sum of connected and disconnected parts: 






DC 



(E„ - Et) («“i‘ - (“t‘) , 




where the last term comes from considering 



-l/2\{Eo-Ei) ($g’|Tf|$o) 




in the equation for doubly excited amplitudes. 



APPENDIX B: 

COMPUTATIONAL IMPLEMENTATION 
Using an extant linear polarizability code, it is straightforward to modify it in 
order to calculate hyperpolarizabilities. The only problem lies in evaluation of the 
perturbation transition matrix between two excited states; 



(hki|h) = 



(s|gi|s)(s|9i|d) 
{d\qi\s) {d\qi\d) 



( 11 ) 



Above matrix is similar to the (h|i/|h) matrix used widely in the EOM theory. 
Elements of this matrix are evaluated as follows: 



(S|gA|S) =<M(ij)-4M(a6), 



(12) 



(S|%|D) = 






(D|9i|S) -t^M {ij,kb} - tfM {ab,cj) , 



( 13 ) 



( 14 ) 



(D\qt\D) = M Uk) - a M (be) . 



u 






( 15 ) 



The algorithm is constructed in such a way, that the whole transition matrix is 
never explicitly evaluated nor stored in the memory or disk. Instead, only one row 
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of this matrix is being formed, and then it is multiplied by a intermediate vector 
(h|Xjt±|0) to create one element of the next order intermediate; vector (h|Xfc±|0) 
to create one element of the next order intermediate; 

I 

( 16 ) 



where the intermediate X amplitudes are defined as 



(h|X*.±|0) = Ro (ioi) * (h|gfc|0). 




The final intermediate, order of which depends on the order of the nonlinear 
process we want to calculate, is multiplied by the left vector containing lambda: 

xl”’ = ( 0 | (l+A),-|h)»(h|4r'*|0> 




and the process is repeated for all transition moments. 



APPENDIX C: 

NAMES, ABBREVIATIONS AND UNITS 

Table 1: Names and abbreviations for various linear and nonlinear processes 



Property 


Process 


Abbreviation 


/3(0;a;, -w) 


Optical Rectification 


OR 


(5{—2uj\ Lo, cj) 


Second Harmonic Generation 


SHG 




dc-Pockels Effect 


dc-P 


7 (— O’; u), LO, —Lo) 


Intensity Dependent Refractive Index 


IDRI 


7(— o)i; a>i, 0)2) —^ 2 ) 


ac-Kerr Effect 


ac-K 


7 (— 3u;; uo, lo, lo) 


Third Harmonic Generation 


THG 


7(-u;;w,0,0) 


dc-Kerr Effect 


dc-K 


7(0,w, -w,0) 


dc-Optical Rectification 


dc-OR 


7(— 2a>; lo, lo, 0) 


dc-Second Harmonic Generation 


dc-SHG 


i^(—Lo; LO, 0, 0) 


Gotton-Mouton Effect 


GME 



Table 2: Units and Equivalences 



Property 


Atomic Unit 


SI Equivalent 


cgs-esu Equivalent 


a 


^2^2 27— 1 


1.6488 X 10-41 C^m^J-i 


2.9635 X 10“^^ Fr^cm^erg-^ 




^3^3 rp—2 


3.20636 X 10-^^ C^mM-2 


8.63922 X Fr^cm^erg"^ 


7 


^ 4^4 ZTi — 3 


6.23538 X 10-®^ G^m^J-^ 


5.03670 X Fr^cm'^erg~^ 



98 



REFERENCES 



[1] P.A. Franken, A.E. Hill, C.W. Peters, and G. Weinreich. Phys. Rev. Lett., 
7:118, 1961. 

[2] M. Kaiser and C.G.B. Garrett. Phys. Rev. Lett., 7:229, 1961. 

[3] H. Sekino and R.J. Bartlett. J. Chem. Phys., 85:976, 1986. 

[4] H. Sekino and R.J. Bartlett. Int. J. Quantum Chem., 43:119-199, 1992. 

[5] H. Sekino and R.J. Bartlett. J. Chem. Phys., 94(5):3665, 1991. 

[6] P. J0rgensen. Ann. Rev. Phys. Chem., 84:359, 1975. 

[7] V. Spirko, Y. Luo, H. Agren, and P. Jprgensen. J. Chem. Phys., 99:9815- 
9815, 1993. 

[8] D.M. Bishop. Rev. Mod. Phys., 62:343, 1990. 

[9] P. Jprgensen and T. Helgaker. J. Chem. Phys., 89:3654, 1988. 

[10] Y. Luo, H. Agren, P. Jprgensen, and K.V. Mikkelsen. Adv. Quantum Chem., 
26:165, 1995. 

[11] P. Albertsen, P. Jprgensen, and D.L. Yeager. Mol. Phys., 41:409, 1980. 

[12] M. Jaszunski, A. Rizzo, and D.L. Yeager. Chem. Phys. Lett., 149:79, 1988. 

[13] J.E. Rice and N.C. Handy. J. Chem. Phys., 94:4959, 1991. 

[14] J.E. Rice and N. G. Handy. Int. J. Quantum Chem., 43:91, 1992. 

[15] K. Sasagane, F. Aiga, and R. Itoh. J. Chem. Phys., 99:3738, 1993. 

[16] F. Aiga, K. Sasagane, and R. Itoh. J. Chem. Phys., 99:3779, 1993. 

[17] F. Aiga, K. Sasagane, and R. Itoh. Int. J. Quantum Chem., 51:87, 1994. 

[18] K. Emrich. Nucl. Phys., A351(l):379, 1981. 

[19] D.C. Comeau and R.J. Bartlett. Chem. Phys. Lett., 207(4):414, 1993. 

[20] J.F. Stanton and R.J. Bartlett. J. Chem. Phys., 97(4):2299, 1993. 

[21] J. F. Stanton and R. J. Bartlett. J. Chem. Phys., 99:5178, 1993. 



99 



100 



[22] A.D. Buckingham. Adv. Chem. Phys., 12:107, 1967. 

[23] D.A. Kleinman. Phys. Rev. A, 126:1977, 1962. 

[24] T. Inoue and S. Iwata. Chem. Phys. Lett., 167(6):566, 1990. 

[25] T.-S. Nee, R.G. Parr, and R.J. Bartlett. J. Chem. Phys., 64:2216, 1976. 

[26] C.E. Dykstra and P.G. Jasien. Chem. Phys. Lett., 109:388, 1984. 

[27] D. N. Zubarev. Soviet Physics Uspekhi, 3:1870, 1960. 

• • 

[28] J. Linderberg and Y. Ohrn. Propagators in Quantum Chemistry. Academic 
Press, London and New York, 1973. 

[29] P. J0rgensen, H. J. Aa Jensen, and J. Olsen. J. Chem. Phys., 82:3235, 1985. 

[30] J. Oddershede, P. J0rgensen, and D. L. Yeager. Comp. Phys. Rep., 2(1):33, 
1989. 

[31] D.C. Hanna, M.A. Yuratich, and D. Cotter. Nonlinear Optics of Free Atoms 
and Molecules. Springer- Verlag, New York, 1979. 

[32] D.M. Bishop. J. Chem. Phys., 100:6536, 1994. 

[33] D.J. Thouless. The Quantum Mechanics of Many-Body Systems. Academic, 
New York, 1961. 

[34] D. L. Yeager and P. Jprgensen. Chem. Phys. Lett., 65:77, 1979. 

[35] E. Dalgaard. J. Chem. Phys., 72:816, 1980. 

[36] E.S. Nielsen, P. Jprgensen, and J. Oddershede. J. Chem. Phys., 73:6238, 
1980. 

[37] H. J. Monkhorst. Int. J. Quantum Chem. Symp., 11:421, 1977. 

[38] E. Dalgaard and H.J. Monkhorst. Phys. Rev. A, 28:1217, 1983. 

[39] H. Koch and P. Jprgensen. J. Chem. Phys., 93:3333, 1990. 

[40] C. Hattig, O. Christiansen, H. Koch, and P. Jprgensen. Chem. Phys. Lett., 
269:428, 1997. 

[41] I. Shavitt and R.J. Bartlett, in preparation. 

[42] E. A. Salter, Gary W. Trucks, and Rodney J. Bartlett. J. Chem. Phys., 
90:1752, 1989. 

[43] R.J. Bartlett. In P. Jprgensen and J. Simons, editors. Geometrical Derivatives 
of Energy surfaces and Geometrical Properties. Reidel, Dordrecht, Holland, 
1986. 



101 



[44] H. Sekino and R. J. Bartlett. Chem. Phys. Lett., 225:486, 1994. 

[45] R.W. Boyd. Nonlinear Optics. Academic Press, New York, 1992. 

[46] P.B. Rozyczko, P. Piecuch, and R.J. Bartlett, to be published. 

[47] H. Sekino and R. J. Bartlett. Chem. Phys. Lett., 234:87, 1995. 

[48] P. B. Rozyczko and R. J. Bartlett. J. Chem. Phys., 107:10823, 1997. 

[49] L. Brillouin. J. Phys. Radium, 3:373, 1932. 

[50] E. Wigner. Math. Naturw. Anz. Ungar. Akad. Wiss., 53:477, 1935. 

[51] I. Lindgren. J. Phys. B., 7:2441, 1974. 

[52] I. Lindgren and J. Morrison. Atomic Many-Body Theory. Springer, Berlin, 
1982. 

[53] P. B. Rozyczko, S. A. Perera, M. Nooijen, and R. J. Bartlett. J. Chem. 
Phys., 107:6736, 1997. 

[54] P.-O. Lbwdin. Phys. Rev. A, A139:357, 1965. 

[55] P.-O. Lowdin. J. Chem. Phys., 43:175, 1965. 

[56] P.-O. Lowdin. Int. J. Quantum Chem., 21:69, 1981. 

[57] E. Brandas and R.J. Bartlett. Chem. Phys. Lett., 8:153, 1981. 

[58] E. Brandas and R.J. Bartlett. J. Chem. Phys., 56:5467, 1972. 

[59] I. Hubac and P. Neogrady. Phys. Rev. A, A50:4558, 1994. 

[60] P. Carsky, V. Hrouda V., Sychrovsky, I. Hubac, P. Babinec, P. Mach, J. Ur- 
ban, and J. Masik. Collect. Czech. Chem. Commun., 60:1419, 1995. 

[61] J. Masik and I. Hubac. Collect. Czech. Chem. Commun., 62:829, 1997. 

[62] P.L. Knight and L. Allen. Phys. Lett., 38A:99, 1972. 

[63] J. Oddershede. Adv. Quantum Chem., 11(1):275, 1978. 

[64] P. G. Szalay, M. Nooijen, and R. J. Bartlett. J. Chem. Phys., 93:281, 1995. 

[65] L. Meissner and R. J. Bartlett. J. Chem. Phys., 92:561, 1990. 

[66] R. Kobayashi, H. Koch, and Poul Jprgensen. Chem. Phys. Lett, 219(5) :30, 
1994. 

[67] H. Sekino and R.J. Bartlett, to be published. 



102 



[68] J. Geertsen, M. Rittby, and R.J. Bartlett. Chem. Phys. Lett., 164(1):57, 
1989. 

[69] S. R. Gwaltney, M. Nooijen, and R. J. Bartlett. Chem. Phys. Lett., 248:189, 

1996. 

[70] S. A. Perera, M. Nooijen, and R. J. Bartlett. Chem. Phys. Lett., 266:456, 

1997. 

[71] A.J. Sadlej. Collect. Czech. Chem. Commun., 53:1995, 1988. 

[72] A.J. Sadlej. Theor. Chim. Acta, 79:123, 1991. 

[73] H. Sekino and R.J. Bartlett. Chem. Phys. Lett., 234:87, 1995. 

[74] Aces II program is a product of the Quantum Theory Project, University 
of Florida. Authors: J. F. Stanton, J. Gauss, J. D. Watts, M. Nooijen, 
N. Oliphant, , S. A. Perera, P. G. Szalay, W. J. Lauderdale, S. R. Gwaltney, 
S. Beck, A. Baikova, D. E. Bernholdt, K. K. Baeck, H. Sekino, P. Rozyczko, 
C. Huber, R. J. Bartlett. Integral packages included are VMOL (J. Almlof, 
J. Olsen P. Taylor); VPROPS (P. R. Taylor); modified version of ABACUS 
integral derivative package (T. U. Helgaker, H. J. Aa. Jensen, P. Jprgensen, 
and P. R. Taylor.). 

[75] M.A. Spackman. J. Phys. Chem., 93:7594, 1989. 

[76] R. D. Amos, N. C. Handy, P. J. Knowles, J. E. Rice, and A. J. Stone. J. 
Phys. Chem., 89:2186, 1985. 

[77] R.J. Bartlett and J.E. Stanton. In K. B. Lipkowitz and D. B. Boyd, edi- 
tors, Reviews in Computational Chemistry, volume 5, pages 65-169. VCH 
Publishers, Inc., New York, 1994. 

[78] H. Sekino and R. J. Bartlett. J. Chem. Phys., 98:3022, 1993. 

[79] J. Oddershede and E. N. Svendsen. Chem. Phys., 64:359, 1982. 

[80] E.A. Salter, H. Sekino, and R.J. Bartlett. J. Chem. Phys., 87:502, 1987. 

[81] E. W. Washburn. International Critical Tables. McGraw-Hill, New York, 
1982. 

[82] G. R. Alms, A. K. Burnham, and W. H. Elygare. J. Chem. Phys., 63:3321, 
1975. 

[83] N. J. Bridge and A. D. Buckingham. Proc. Roy. Soc., Ser. A, A295:334, 
1966. 

[84] Y. Luo, O. Vahtras, H. Agren, and P. Jprgensen. Chem. Phys. Lett., 205:555, 
1997. 



103 



[85] R.J. Bartlett and H. Sekino. In S.P. Kama and A.T. Yates, editors, Non- 
Linear Optical Materials: Theory and Modelling, page 23. American Chem- 
ical Society, Boca Raton, 1996. 

[86] P. Norman, Y. Luo, D. Jonsson, and H. Agren. J. Chem. Phys., 106:1827, 
1997. 

[87] In M.J. Astle R.C. Weast and W.H. Beyer, editors. Handbook of Chemistry 
and Physics, pages E-70. CRC, Boca Raton, 1985. 

[88] J. F. Stanton. Phys. Rev. A, A49:1698, 1994. 

[89] W. Rijks and P. E. S. Wormer. J. Chem. Phys., 88(9):5704, 1988. 

[90] W. Rijks and P. E. S. Wormer. J. Chem. Phys., 90(11):6507, 1989. 

[91] M. A. Spackman. J. Chem. Phys., 94:1295, 1991. 

[92] D. J. Margoliash and W. J. Meath. J. Chem. Phys., 68:1426, 1978. 

[93] G. Starkschall and R. G. Gordon. J. Chem. Phys., 54:6635, 1971. 

[94] P. W. Fowler, P. Jorgensen, and J. Olsen. J. Chem. Phys., 93:7256, 1990. 

[95] J. D. Watts and R. J. Bartlett. Chem. Phys. Lett., 258:581, 1996. 

[96] H. Sekino and Rodney J. Bartlett. J. Chem. Phys., 98:3022, 1993. 

[97] R. J. Bartlett and G. D. Purvis III. Phys. Rev. A, 20:1313, 1979. 

[98] H. Sekino and R.J. Bartlett. J. Chem. Phys., 84:2726, 1986. 

[99] D.P. Chong and S.R. Langhoff. J. Chem. Phys., 93:570, 1990. 

[100] M. Jaszunski, P. Jorgensen, and H.J.Aa. Jensen. Chem. Phys. Lett., 191:293, 
1992. 

[101] H. Sekino and R. J. Bartlett. J. Chem. Phys., 98:3022, 1993. 

[102] D.P. Shelton and J.E. Rice. Chem. Rev., 29:3, 1993. 

[103] M.G. Papadopoulos, J. Waite, and A.D. Buckingham. J. Chem. Phys., 
102:371, 1995. 

[104] J.W. Dudley II and J.F. Ward. J. Chem. Phys., 82:4673, 1985. 

[105] D.E. Woon and T.H. Dunning. J. Chem. Phys., 100:2975, 1994. 

[106] S. A. Perera, M. Nooijen, and R. J. Bartlett. J. Chem. Phys., 104:3290, 
1996. 

[107] D.M. Bishop and B. Kirtman. J. Chem. Phys., 95:2646, 1991. 



104 



[108] S.M. Cybulski and D.M. Bishop. J. Chem. Phys., 100:2019, 1994. 

[109] C. Hattig and B.A. Hess. Chem. Phys. Lett., 233:359, 1995. 

[110] T. C. Caves and M. Karpins. J. Chem. Phys., 50:3649, 1969. 

[111] R.J. Bartlett. In D. R. Yarkony, editor. Modern Electronic Structure Theory, 
part I, page 1047. CRC, Boca Raton, 1995. 

[112] P. Norman, Y. Luo, D. Jonsson, and H. Agren. J. Chem. Phys., 106:1827, 
1997. 

[113] G. J. B. Hurst, M. Dupuis, and E. Clementi. J. Chem. Phys., 89:385, 1988. 

[114] W. Haugen and M. Traetteberg. Acta Chem. Scand., 20:1726, 1966. 

[115] D. M. Bishop. Chem. Phys. Lett., 69:5438, 1978. 

[116] D. P. Shelton and J. E. Rice. Chem. Rev., 29:3, 1993. 

[117] J. E. Ward and D. S. Elliott. J. Chem. Phys., 69:5438, 1978. 

[118] B. Kirtman, J. T. Toto, C. Breneman, C. P. de Melo, and D. M. Bishop. J. 
Chem. Phys., in press. 

[119] J. D. Watts, S. R. Gwaltney, and R. J. Bartlett. J. Chem. Phys., 105:6979, 
1996. 

[120] J. D. Watts. Private Communication. 

[121] R. McDiarmid. J. Chem. Phys., 64:514, 1976. 

[122] L. Serrano-Andres, M. Merchan, I. Nebot-Gil, R. Lindh, and B. O. Roos. J. 
Chem. Phys., 98:3151, 1993. 

[123] P. G. Szalay and R. J. Bartlett. J. Chem. Phys., 101:4936, 1994. 

[124] K. Wilson, P. Piecuch, S. Kucharski, and R.J. Bartlett. J. Chem. Phys., to 
be published. 

[125] C. Rizzo, A. Rizzo, and D.M. Bishop. Inti Rev. Phys. Chem., 16:81, 1997. 

[126] S. Coriani, K. Ruud A. Rizzo, and T. Helgaker. Chem. Phys., 216:53, 1997. 

[127] S.A. Perera, H. Sekino, and R.J. Bartlett. J. Chem. Phys., 101:2186, 1994. 

[128] S.A. Perera, M. Nooijen, and R.J. Bartlett. J. Chem. Phys., 104:3290, 1996. 

[129] J.Q. Sun and R.J. Bartlett. J. Chem. Phys., 107:5058, 1997. 

[130] Q.L. Zhou, J.R. Heflin, K.Y. Wong, O. Zamani-Khamiri, and A.F. Garito. 
Phys. Rev. A, 43:1673, 1991. 



105 



[131] D.C. Rodenberger, J.R. Heflin, and A.F. Garito. Phys. Rev. A, 51:3234, 

1995. 

[132] D. Jonsson, Y. Luo P. Norman, and H. Agren. J. Chem. Phys., 105:581, 

1996. 



BIOGRAPHICAL SKETCH 

Piotr Rozyczko was born on December 18, 1968, in a really small city in the 
southeastern part of Poland. He earned his Technician in Analytical Chemistry 
diploma in 1988, being 18. Next step in his career was the Master of Science degree 
in 1993 at the Wroclaw University, where he successfully defended his M.S. thesis 
on the electronic properties of phosphorus analogs of several nitrogen containing 
molecules. 

In August, 1993, he arrived in Gainesville and begun graduate studies at the 
University of Florida. His main interest has been advanced quantum chemistry 
methods, and Quantum Theory Project and especially Dr. Rodney Bartlett’s 
group suited his expectations best. 

Assuming that this thesis is turned in on time, he will graduate in December 
of 1998. 



106 



I certify that I have read this study and that in my opinion it conforms to 
acceptable standards of scholarly presentation and is fully adequate, in scope and 
quality, as a dissertation for the degree of Doctor of Philosophy. 




Graduate Research Professor of 
Chemistry 



I certify that I have read this study and that in my opinion it conforms to 
acceptable standards of scholarly presentation and is fully adequate, in scope and 
quality, as a dissertation for the degree of Doctor of Philosophy. 




N. Yngve Of 
Professor of 




hemistry and Physics 



I certify that I have read this study and that in my opinion it conforms to 
acceptable standards of scholarly presentation and is fully adequate, in scope and 
quality, as a dissertation for the degree of Doctor of Philosophy. 




Professor of Chemistry and Physics 



I certify that I have read this study and that in my opinion it conforms to 



acceptable standards of scholarly presentation 
quality, as a dissertation for the degree o 




is fully adequate, in scope and 
(jf Philosophy. 



HendriK J. Monkhorst 
Professor of Physics 



I certify that I have read this study and that in my opinion it conforms to 
acceptable standards of scholarly presentation and is fully adequate, in scope and 
quality, as a dissertation for the degree of Doctor of Philosophy. 





fjAmes Deyrup 
Professor of Chemistry 



This dissertation was submitted to the Graduate Faculty of the Department 
of Chemistry in the College of Liberal Arts and Sciences, and to the Graduate 
School and was accepted as partial fulfillment of the requirements for the degree 
of Doctor of Philosophy. 

December 1998 

Dean, Graduate School 



